/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/
/*
 */

#include <dxconfig.h>

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <dx/dx.h>
#include "_connectgrids.h"

struct arg_radius
{
  Object ino;
  Field base;
  float radius;
  float exponent;
  Array missing;
};

struct arg_nearest
{
  Object ino;
  Field base;
  int number;
  float *radius;
  float radiusvalue;
  float exponent;
  Array missing;
};

struct arg_scatter
{
  Object ino;
  Field base;
  Array missing;
};

typedef struct matrix2D
{
  float A[2][2];
  float b[2];
} Matrix2D;

typedef struct vector2D
{
  float x, y;
} Vector2D;

#define PI 3.14156
#define E .0001 /* fuzz */
#define FLOOR( x )                       \
  ( ( x ) + E > 0 ? ( int )( ( x ) + E ) \
                  : ( int )( ( x ) + E ) - 1 ) /* fuzzy floor */
#define CEIL( x )                            \
  ( ( x ) - E > 0 ? ( int )( ( x ) - E ) + 1 \
                  : ( int )( ( x ) - E ) ) /* fuzzy ceiling */
#define ABS( x ) ( ( x ) < 0 ? -( x ) : ( x ) )
#define RND( x ) ( x < 0 ? x - 0.5 : x + 0.5 )
#define POW( x, y ) ( pow( (double)x, (double)y ) )

#define INSERT_NEAREST_BUFFER( TYPE )                                         \
  {                                                                           \
    TYPE *holder = (TYPE *)bufdata;                                           \
    float distancesq = 0;                                                     \
    switch ( shape_base[0] )                                                  \
    {                                                                         \
      case( 1 ) :                                                             \
        distancesq = ( gridposition1D - inputposition_ptr[l] ) *              \
                     ( gridposition1D - inputposition_ptr[l] );               \
        break;                                                                \
      case( 2 ) :                                                             \
        distancesq = ( gridposition2D.x - inputposition_ptr[2 * l] ) *        \
                         ( gridposition2D.x - inputposition_ptr[2 * l] ) +    \
                     ( gridposition2D.y - inputposition_ptr[2 * l + 1] ) *    \
                         ( gridposition2D.y - inputposition_ptr[2 * l + 1] ); \
        break;                                                                \
      case( 3 ) :                                                             \
        distancesq = ( gridposition.x - inputposition_ptr[3 * l] ) *          \
                         ( gridposition.x - inputposition_ptr[3 * l] ) +      \
                     ( gridposition.y - inputposition_ptr[3 * l + 1] ) *      \
                         ( gridposition.y - inputposition_ptr[3 * l + 1] ) +  \
                     ( gridposition.z - inputposition_ptr[3 * l + 2] ) *      \
                         ( gridposition.z - inputposition_ptr[3 * l + 2] );   \
        break;                                                                \
    }                                                                         \
    /* first try to insert into the already present buffer */                 \
    for ( m = 0; m < numvalid; m++ )                                          \
    {                                                                         \
      if ( distancesq < bufdistance[m] )                                      \
      {                                                                       \
        if ( numvalid < numnearest )                                          \
          numvalid++;                                                         \
        for ( n = numvalid - 1; n > m; n-- )                                  \
        {                                                                     \
          bufdistance[n] = bufdistance[n - 1];                                \
          memcpy( &holder[shape[0] * n], &holder[shape[0] * ( n - 1 )],       \
                  DXGetItemSize( newcomponent ) );                            \
        }                                                                     \
        bufdistance[m] = distancesq;                                          \
        memcpy( &holder[shape[0] * m], tmpdata,                               \
                DXGetItemSize( newcomponent ) );                              \
        goto filled_buf##TYPE;                                                \
      }                                                                       \
    }                                                                         \
    if ( numvalid == 0 )                                                      \
    {                                                                         \
      bufdistance[0] = distancesq;                                            \
      memcpy( &holder[0], tmpdata, DXGetItemSize( newcomponent ) );           \
      numvalid++;                                                             \
    }                                                                         \
    filled_buf##TYPE : tmpdata += datastep;                                   \
    continue;                                                                 \
  }

#define ACCUMULATE_DATA( TYPE )                                         \
  {                                                                     \
    TYPE *holder;                                                       \
    float dis;                                                          \
    float radiussq = radiusvalue * radiusvalue;                         \
    int dc;                                                             \
    for ( ii = 0, holder = (TYPE *)bufdata; ii < numvalid; ii++ )       \
    {                                                                   \
      dis = sqrt( bufdistance[ii] );                                    \
      if ( bufdistance[ii] != 0.0 )                                     \
      {                                                                 \
        if ( radius )                                                   \
        {                                                               \
          if ( bufdistance[ii] <= radiussq )                            \
          {                                                             \
            distanceweight += 1.0 / POW( dis, exponent );               \
            for ( dc = 0; dc < shape[0]; dc++ )                         \
            {                                                           \
              data[dc] += (float)( holder[dc] ) / POW( dis, exponent ); \
            }                                                           \
          }                                                             \
        }                                                               \
        else                                                            \
        {                                                               \
          distanceweight += 1.0 / ( POW( dis, exponent ) );             \
          for ( dc = 0; dc < shape[0]; dc++ )                           \
          {                                                             \
            data[dc] += (float)holder[dc] / POW( dis, exponent );       \
          }                                                             \
        }                                                               \
      }                                                                 \
      else                                                              \
      {                                                                 \
        distanceweight = 1.0;                                           \
        for ( dc = 0; dc < shape[0]; dc++ )                             \
        {                                                               \
          data[dc] = (float)holder[dc];                                 \
        }                                                               \
        break;                                                          \
      }                                                                 \
      holder += shape[0];                                               \
    }                                                                   \
  }

#define CHECK_AND_ALLOCATE_MISSING( missing, size, type, type_enum, shape,     \
                                    ptr )                                      \
  {                                                                            \
    ptr = (type *)DXAllocateLocal( size );                                     \
    if ( !ptr )                                                                \
      goto error;                                                              \
    if ( missing )                                                             \
    {                                                                          \
      if ( !DXExtractParameter( (Object)missing, type_enum, shape, 1,          \
                                (Pointer)ptr ) )                               \
      {                                                                        \
        DXSetError( ERROR_DATA_INVALID, "missing must match all "              \
                                        "position-dependent components; does " \
                                        "not match %s",                        \
                    name );                                                    \
        goto error;                                                            \
      }                                                                        \
    }                                                                          \
    else                                                                       \
    {                                                                          \
      DXWarning( "Grid missing values assigned to 0 and marked invalid" );     \
      for ( i = 0; i < shape; i++ )                                            \
        ptr[i] = 0;                                                            \
    }                                                                          \
  }

static Error ConnectNearest( Group, Field, int, float *, float, float, Array );
static Error FillBuffersIrreg( int, int, ubyte *, float *, float *, float *,
                               InvalidComponentHandle, InvalidComponentHandle,
                               int *, int, float );
static Error FillDataArray( Type, Array, Array, ubyte *, float *, int, int, int,
                            float, Array, char * );
static Error ConnectRadiusField( Pointer );
static Error ConnectNearestField( Pointer );
static Error ConnectScatterField( Pointer );
static Error ConnectScatter( Field, Field, Array );
static Error ConnectRadius( Field, Field, float, float, Array );
static Error ConnectRadiusRegular( Field, Field, float, float, Array );
static Error ConnectRadiusIrregular( Field, Field, float, float, Array );

static Vector2D Vec2D( double x, double y )
{
  Vector2D v;
  v.x = x;
  v.y = y;
  return v;
}

static Vector2D Apply2D( Vector2D vec, Matrix2D mat )
{
  Vector2D result;

  result.x = vec.x * mat.A[0][0] + vec.y * mat.A[1][0] + mat.b[0];
  result.y = vec.x * mat.A[0][1] + vec.y * mat.A[1][1] + mat.b[1];
  return result;
}

static Matrix2D Invert2D( Matrix2D matrix2d )
{
  float a, b, c, d, notsing;
  Matrix2D invert;
  Vector2D vec2d, invertb;

  a = matrix2d.A[0][0];
  b = matrix2d.A[0][1];
  c = matrix2d.A[1][0];
  d = matrix2d.A[1][1];

  /* check singularity */
  notsing = a * d - b * c;
  if ( notsing )
  {
    invert.A[0][0] = ( 1.0 / notsing ) * ( d );
    invert.A[0][1] = ( 1.0 / notsing ) * ( -b );
    invert.A[1][0] = ( 1.0 / notsing ) * ( -c );
    invert.A[1][1] = ( 1.0 / notsing ) * ( a );
    invert.b[0] = 0;
    invert.b[1] = 0;
  }

  invertb.x = matrix2d.b[0];
  invertb.y = matrix2d.b[1];
  vec2d = Apply2D( invertb, invert );
  invert.b[0] = -vec2d.x;
  invert.b[1] = -vec2d.y;
  return invert;
}

Error _dxfConnectNearestObject( Object ino, Object base, int numnearest,
                                float *radius, float exponent, Array missing )
{
  /* recurse on base to the field level */
  struct arg_nearest arg;
  Object subbase;
  int i;

  switch ( DXGetObjectClass( base ) )
  {
    case CLASS_FIELD:
      /* create the task group */
      arg.ino = ino;
      arg.base = (Field)base;
      arg.number = numnearest;
      arg.radius = radius;
      arg.exponent = exponent;
      if ( radius )
        arg.radiusvalue = *radius;
      else
        arg.radiusvalue = 0.0;
      arg.missing = missing;
      if ( !DXAddTask( ConnectNearestField, ( Pointer ) & arg, sizeof( arg ),
                       0.0 ) )
        goto error;
      break;
    case CLASS_GROUP:
      for ( i = 0; ( subbase = DXGetEnumeratedMember( (Group)base, i, NULL ) );
            i++ )
      {
        if ( !( _dxfConnectNearestObject( (Object)ino, (Object)subbase,
                                          numnearest, radius, exponent,
                                          missing ) ) )
          goto error;
      }
      break;
    default:
      DXSetError( ERROR_DATA_INVALID, "base must be a field or a group" );
      goto error;
  }
  return OK;

error:
  return ERROR;
}

Error _dxfConnectRadiusObject( Object ino, Object base, float radius,
                               float exponent, Array missing )
/* recurse on the base to the field level */
{
  struct arg_radius arg;
  Object subbase;
  int i;

  switch ( DXGetObjectClass( base ) )
  {
    case CLASS_FIELD:
      /* create the task group */
      arg.ino = ino;
      arg.base = (Field)base;
      arg.radius = radius;
      arg.exponent = exponent;
      arg.missing = missing;
      if ( !DXAddTask( ConnectRadiusField, ( Pointer ) & arg, sizeof( arg ),
                       0.0 ) )
        goto error;
      break;
    case CLASS_GROUP:
      for ( i = 0; ( subbase = DXGetEnumeratedMember( (Group)base, i, NULL ) );
            i++ )
      {
        if ( !( _dxfConnectRadiusObject( (Object)ino, (Object)subbase, radius,
                                         exponent, missing ) ) )
          goto error;
      }
      break;
    default:
      DXSetError( ERROR_DATA_INVALID, "base must be a field or a group" );
      goto error;
  }
  return OK;

error:
  return ERROR;
}

static Error ConnectNearestField( Pointer p )
{
  struct arg_nearest *arg = (struct arg_nearest *)p;
  Object ino, subino, newino = NULL;
  int i, numnearest;
  Field base;
  float *radius;
  float radiusvalue;
  float exponent;
  Array missing;

  ino = arg->ino;
  newino = DXApplyTransform( ino, NULL );

  switch ( DXGetObjectClass( newino ) )
  {

    case CLASS_FIELD:

      base = arg->base;
      radius = arg->radius;
      radiusvalue = arg->radiusvalue;
      missing = arg->missing;
      numnearest = arg->number;
      exponent = arg->exponent;
      if ( !ConnectNearest( (Group)newino, base, numnearest, radius,
                            radiusvalue, exponent, missing ) )
        goto error;
      break;

    case CLASS_GROUP:

      switch ( DXGetGroupClass( (Group)newino ) )
      {
        case CLASS_COMPOSITEFIELD:
          /* for now, can't handle a composite field */
          DXSetError( ERROR_NOT_IMPLEMENTED,
                      "input cannot be a composite field" );
          goto error;
          base = arg->base;
          radius = arg->radius;
          radiusvalue = arg->radiusvalue;
          missing = arg->missing;
          numnearest = arg->number;
          exponent = arg->exponent;
          if ( !ConnectNearest( (Group)newino, base, numnearest, radius,
                                radiusvalue, exponent, missing ) )
            goto error;
          break;
        default:
          DXSetError( ERROR_NOT_IMPLEMENTED, "input must be a single field" );
          goto error;
          /* recurse on the members of the group */
          for ( i = 0;
                ( subino = DXGetEnumeratedMember( (Group)newino, i, NULL ) );
                i++ )
          {
            base = arg->base;
            radius = arg->radius;
            radiusvalue = arg->radiusvalue;
            missing = arg->missing;
            numnearest = arg->number;
            exponent = arg->exponent;
            if ( !ConnectNearest( (Group)subino, base, numnearest, radius,
                                  radiusvalue, exponent, missing ) )
              goto error;
          }
          break;
      }
      break;
    default:
      DXSetError( ERROR_DATA_INVALID, "input must be a field or a group" );
      goto error;
  }
  DXDelete( (Object)newino );
  return OK;

error:
  DXDelete( (Object)newino );
  return ERROR;
}

static Error ConnectRadiusField( Pointer p )
{
  struct arg_radius *arg = (struct arg_radius *)p;
  Object ino, newino = NULL;
  Field base;
  float radius;
  float exponent;
  Array missing;

  /* this one recurses on object ino */

  ino = arg->ino;
  newino = DXApplyTransform( ino, NULL );

  switch ( DXGetObjectClass( newino ) )
  {

    case CLASS_FIELD:
      base = arg->base;
      radius = arg->radius;
      missing = arg->missing;
      exponent = arg->exponent;
      if ( !ConnectRadius( (Field)newino, base, radius, exponent, missing ) )
        goto error;
      break;

    case CLASS_GROUP:
      DXSetError( ERROR_NOT_IMPLEMENTED, "input cannot be a group" );
      goto error;
      break;

    default:
      DXSetError( ERROR_DATA_INVALID, "input must be a field" );
      goto error;
  }
  DXDelete( (Object)newino );
  return OK;

error:
  DXDelete( (Object)newino );
  return ERROR;
}

static Error ConnectNearest( Group ino, Field base, int numnearest,
                             float *radius, float radiusvalue, float exponent,
                             Array missing )
{
  int componentsdone, anyinvalid = 0;
  Matrix matrix;
  Matrix2D matrix2D;
  Array positions_base, positions_ino, oldcomponent = NULL, newcomponent;
  InvalidComponentHandle out_invalidhandle = NULL, in_invalidhandle = NULL;
  int counts[8], rank_base, shape_base[8], shape_ino[8], rank_ino;
  int rank, shape[8], i, j, k, l, m, n, numvalid, pos_count = 0;
  int numitemsbase, compcount, numitemsino, numitems;
  float distanceweight, distance = 0, *data = NULL, *dp, *scratch = NULL;
  ArrayHandle handle = NULL;
  Type type;
  Point gridposition = {0, 0, 0};
  Vector2D gridposition2D = {0, 0};
  float gridposition1D = 0, origin = 0, delta = 0;
  float *bufdistance = NULL;
  char *bufdata = NULL;
  char *name;
  Category category;
  float *inputposition_ptr;
  int data_count, invalid_count = 0, ii, datastep, dc, regular;
  char *old_data_ptr = NULL, *tmpdata = NULL;
  float *new_data_ptr_f = NULL, *missingval_f = NULL;
  int *new_data_ptr_i = NULL, *missingval_i = NULL;
  short *new_data_ptr_s = NULL, *missingval_s = NULL;
  double *new_data_ptr_d = NULL, *missingval_d = NULL;
  byte *new_data_ptr_b = NULL, *missingval_b = NULL;
  ubyte *new_data_ptr_ub = NULL, *missingval_ub = NULL;
  ushort *new_data_ptr_us = NULL, *missingval_us = NULL;
  uint *new_data_ptr_ui = NULL, *missingval_ui = NULL;

  if ( DXEmptyField( (Field)ino ) )
    return OK;

  positions_base = (Array)DXGetComponentValue( (Field)base, "positions" );
  if ( !positions_base )
  {
    DXSetError( ERROR_MISSING_DATA, "base field has no positions" );
    goto error;
  }
  positions_ino = (Array)DXGetComponentValue( (Field)ino, "positions" );
  if ( !positions_ino )
  {
    DXSetError( ERROR_MISSING_DATA, "input field has no positions" );
    goto error;
  }

  DXGetArrayInfo( positions_base, &numitemsbase, &type, &category, &rank_base,
                  shape_base );

  regular = 0;
  if ( ( DXQueryGridPositions( (Array)positions_base, NULL, counts,
                               (float *)&matrix.b, (float *)&matrix.A ) ) )
  {
    regular = 1;
    if ( shape_base[0] == 2 )
    {
      matrix2D.A[0][0] = matrix.A[0][0];
      matrix2D.A[1][0] = matrix.A[0][1];
      matrix2D.A[0][1] = matrix.A[0][2];
      matrix2D.A[1][1] = matrix.A[1][0];
      matrix2D.b[0] = matrix.b[0];
      matrix2D.b[1] = matrix.b[1];
    }
    else if ( shape_base[0] == 1 )
    {
      origin = matrix.b[0];
      delta = matrix.A[0][0];
    }
  }
  else
  {
    /* irregular positions */
    if ( !( handle = DXCreateArrayHandle( positions_base ) ) )
      goto error;
    scratch = DXAllocateLocal( DXGetItemSize( positions_base ) );
    if ( !scratch )
      goto error;
  }

  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "positions component must be type float, category real" );
    goto error;
  }

  out_invalidhandle =
      DXCreateInvalidComponentHandle( (Object)base, NULL, "positions" );
  if ( !out_invalidhandle )
    goto error;
  DXSetAllValid( out_invalidhandle );
  in_invalidhandle =
      DXCreateInvalidComponentHandle( (Object)ino, NULL, "positions" );
  if ( !in_invalidhandle )
    return ERROR;

  DXGetArrayInfo( positions_ino, &numitemsino, &type, &category, &rank_ino,
                  shape_ino );

  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "positions component must be type float, category real" );
    goto error;
  }
  if ( rank_base == 0 )
    shape_base[0] = 1;
  if ( rank_base > 1 )
  {
    DXSetError( ERROR_DATA_INVALID,
                "rank of base positions cannot be greater than 1" );
    goto error;
  }
  if ( rank_ino == 0 )
    shape_ino[0] = 1;
  if ( rank_ino > 1 )
  {
    DXSetError( ERROR_DATA_INVALID,
                "rank of input positions cannot be greater than 1" );
    goto error;
  }

  if ( ( shape_base[0] < 0 ) || ( shape_base[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, and 3D positions for base supported" );
    goto error;
  }

  if ( shape_base[0] != shape_ino[0] )
  {
    DXSetError( ERROR_DATA_INVALID, "dimensionality of base (%dD) does not "
                                    "match dimensionality of input (%dD)",
                shape_base[0], shape_ino[0] );
    goto error;
  }

  if ( ( shape_ino[0] < 0 ) || ( shape_ino[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, and 3D positions for input supported" );
    goto error;
  }
  /* first set the extra counts to one to simplify the code for handling
     1, 2, and 3 dimensions */
  if ( regular )
  {
    if ( shape_base[0] == 1 )
    {
      counts[1] = 1;
      counts[2] = 1;
    }
    else if ( shape_base[0] == 2 )
    {
      counts[2] = 1;
    }
  }
  else
  {
    counts[0] = numitemsbase;
    counts[1] = 1;
    counts[2] = 1;
  }

  /* need to get the components we are going to copy into the new field */
  /* copied wholesale from gda */
  compcount = 0;
  componentsdone = 0;
  while ( NULL != ( oldcomponent = (Array)DXGetEnumeratedComponentValue(
                        (Field)ino, compcount, &name ) ) )
  {
    Object attr;

    data_count = 0;
    invalid_count = 0;

    if ( DXGetObjectClass( (Object)oldcomponent ) != CLASS_ARRAY )
      goto component_done;

    if ( !strcmp( name, "positions" ) || !strcmp( name, "connections" ) ||
         !strcmp( name, "invalid positions" ) )
      goto component_done;

    /* if the component refs anything, leave it */
    if ( DXGetComponentAttribute( (Field)ino, name, "ref" ) )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "der" );
    if ( attr )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "dep" );
    if ( !attr )
    {
      /* does it der? if so, throw it away */
      attr = DXGetComponentAttribute( (Field)ino, name, "der" );
      if ( attr )
        goto component_done;
      /* it doesn't ref, doesn't dep and doesn't der.
         Copy it to the output and quit */
      if ( !DXSetComponentValue( (Field)base, name, (Object)oldcomponent ) )
        goto error;
      goto component_done;
    }

    if ( DXGetObjectClass( attr ) != CLASS_STRING )
    {
      DXSetError( ERROR_DATA_INVALID, "dependency attribute not type string" );
      goto error;
    }

    if ( strcmp( DXGetString( (String)attr ), "positions" ) )
    {
      DXWarning( "Component '%s' is not dependent on positions!", name );
      DXWarning( "Regrid will remove '%s' and output the base grid 'data' "
                 "component if it exists",
                 name );
      goto component_done;
    }

    /* if the component happens to be "invalid positions" ignore it */
    if ( !strcmp( name, "invalid positions" ) )
      goto component_done;

    componentsdone++;

    DXGetArrayInfo( (Array)oldcomponent, &numitems, &type, &category, &rank,
                    shape );
    if ( rank == 0 )
      shape[0] = 1;

    /* check that the missing component, if present, matches the component
       we're working on */
    if ( missing )
      switch ( type )
      {
        case( TYPE_FLOAT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      float, type, shape[0], missingval_f );
          break;
        case( TYPE_INT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      int, type, shape[0], missingval_i );
          break;
        case( TYPE_DOUBLE ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      double, type, shape[0], missingval_d );
          break;
        case( TYPE_SHORT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      short, type, shape[0], missingval_s );
          break;
        case( TYPE_BYTE ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      byte, type, shape[0], missingval_b );
          break;
        case( TYPE_UINT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      uint, type, shape[0], missingval_ui );
          break;
        case( TYPE_USHORT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      ushort, type, shape[0], missingval_us );
          break;
        case( TYPE_UBYTE ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      ubyte, type, shape[0], missingval_ub );
          break;
        default:
          DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
          goto error;
      }

    if ( rank == 0 )
      shape[0] = 1;
    if ( rank > 1 )
    {
      DXSetError( ERROR_NOT_IMPLEMENTED, "rank larger than 1 not implemented" );
      goto error;
    }
    if ( numitems != numitemsino )
    {
      DXSetError(
          ERROR_BAD_PARAMETER,
          "number of items in %s component not equal to number of positions",
          name );
      goto error;
    }
    newcomponent = DXNewArrayV( type, category, rank, shape );
    if ( !DXAddArrayData( newcomponent, 0, numitemsbase, NULL ) )
      goto error;
    old_data_ptr = (char *)DXGetArrayData( oldcomponent );

    switch ( type )
    {
      case( TYPE_FLOAT ) :
        new_data_ptr_f = (float *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_INT ) :
        new_data_ptr_i = (int *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_SHORT ) :
        new_data_ptr_s = (short *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_BYTE ) :
        new_data_ptr_b = (byte *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_DOUBLE ) :
        new_data_ptr_d = (double *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_UINT ) :
        new_data_ptr_ui = (uint *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_UBYTE ) :
        new_data_ptr_ub = (ubyte *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_USHORT ) :
        new_data_ptr_us = (ushort *)DXGetArrayData( newcomponent );
        break;
      default:
        DXSetError( ERROR_DATA_INVALID, "unrecognized data type" );
        goto error;
    }

    /* need to allocate a buffer of length n to hold the nearest neighbors*/
    bufdata = (char *)DXAllocateLocal(
        numnearest * DXGetItemSize( newcomponent ) * shape[0] );
    if ( !bufdata )
      goto error;
    bufdistance = (float *)DXAllocateLocal( numnearest * sizeof( float ) );
    if ( !bufdistance )
      goto error;
    data = (float *)DXAllocateLocalZero( sizeof( float ) * shape[0] );
    if ( !data )
      goto error;
    datastep = DXGetItemSize( newcomponent );

    /* loop over all the grid points */
    inputposition_ptr = (float *)DXGetArrayData( positions_ino );
    for ( i = 0; i < counts[0]; i++ )
    {
      for ( j = 0; j < counts[1]; j++ )
      {
        for ( k = 0; k < counts[2]; k++ )
        {
          if ( regular )
          {
            switch ( shape_base[0] )
            {
              case( 3 ) :
                gridposition = DXApply( DXVec( i, j, k ), matrix );
                break;
              case( 2 ) :
                gridposition2D = Apply2D( Vec2D( i, j ), matrix2D );
                break;
              case( 1 ) :
                gridposition1D = origin + i * delta;
            }
          }
          else
          {
            switch ( shape_base[0] )
            {
              case( 3 ) :
                dp = DXGetArrayEntry( handle, pos_count, scratch );
                gridposition = DXPt( dp[0], dp[1], dp[2] );
                break;
              case( 2 ) :
                dp = DXGetArrayEntry( handle, pos_count, scratch );
                gridposition2D.x = dp[0];
                gridposition2D.y = dp[1];
                break;
              case( 1 ) :
                dp = DXGetArrayEntry( handle, pos_count, scratch );
                gridposition1D = dp[0];
                break;
            }
            pos_count++;
          }
          for ( ii = 0; ii < numnearest; ii++ )
          {
            bufdistance[ii] = DXD_MAX_FLOAT;
          }
          /* numvalid is how many valid data points are in the buffer
             for this particular grid point */
          numvalid = 0;
          for ( l = 0, tmpdata = old_data_ptr; l < numitemsino; l++ )
          {
            /* don't do anything if the input is invalid */
            if ( DXIsElementValid( in_invalidhandle, l ) )
            {
              switch ( type )
              {
                case( TYPE_FLOAT ) :
                  INSERT_NEAREST_BUFFER( float );
                  break;
                case( TYPE_INT ) :
                  INSERT_NEAREST_BUFFER( int );
                  break;
                case( TYPE_SHORT ) :
                  INSERT_NEAREST_BUFFER( short );
                  break;
                case( TYPE_DOUBLE ) :
                  INSERT_NEAREST_BUFFER( double );
                  break;
                case( TYPE_BYTE ) :
                  INSERT_NEAREST_BUFFER( byte );
                  break;
                case( TYPE_UBYTE ) :
                  INSERT_NEAREST_BUFFER( ubyte );
                  break;
                case( TYPE_USHORT ) :
                  INSERT_NEAREST_BUFFER( ushort );
                  break;
                case( TYPE_UINT ) :
                  INSERT_NEAREST_BUFFER( uint );
                  break;
                default:
                  DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
                  goto error;
              }
            }
          }
          /* have now looped over all the scattered points */
          distanceweight = 0;
          memset( data, 0, shape[0] * sizeof( float ) );

          switch ( type )
          {
            case( TYPE_FLOAT ) :
              ACCUMULATE_DATA( float );
              break;
            case( TYPE_INT ) :
              ACCUMULATE_DATA( int );
              break;
            case( TYPE_BYTE ) :
              ACCUMULATE_DATA( byte );
              break;
            case( TYPE_DOUBLE ) :
              ACCUMULATE_DATA( double );
              break;
            case( TYPE_SHORT ) :
              ACCUMULATE_DATA( short );
              break;
            case( TYPE_UINT ) :
              ACCUMULATE_DATA( uint );
              break;
            case( TYPE_USHORT ) :
              ACCUMULATE_DATA( ushort );
              break;
            case( TYPE_UBYTE ) :
              ACCUMULATE_DATA( ubyte );
              break;
            default:
              DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
              goto error;
          }

          /* now stick the data value into the new data component */
          switch ( type )
          {
            case( TYPE_FLOAT ) :
              if ( distanceweight != 0. )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_f[data_count++] =
                      (float)( data[dc] / distanceweight );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_f[data_count++] = missingval_f[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_f[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_INT ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_i[data_count++] =
                      (int)( RND( data[dc] / distanceweight ) );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_i[data_count++] = missingval_i[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_i[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_SHORT ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_s[data_count++] =
                      (short)( RND( data[dc] / distanceweight ) );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_s[data_count++] = missingval_s[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_s[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_BYTE ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_b[data_count++] =
                      ( byte )( RND( data[dc] / distanceweight ) );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_b[data_count++] = missingval_b[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_b[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_DOUBLE ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_d[data_count++] =
                      (double)( data[dc] / distanceweight );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_d[data_count++] = missingval_d[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_d[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_UBYTE ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_ub[data_count++] =
                      ( ubyte )( RND( data[dc] / distanceweight ) );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_ub[data_count++] = missingval_ub[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_ub[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_USHORT ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_us[data_count++] =
                      ( ushort )( RND( data[dc] / distanceweight ) );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_us[data_count++] = missingval_us[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_us[data_count++] = 0;
                  }
                }
              }
              break;
            case( TYPE_UINT ) :
              if ( distanceweight != 0 )
              {
                for ( dc = 0; dc < shape[0]; dc++ )
                {
                  new_data_ptr_ui[data_count++] =
                      ( uint )( RND( data[dc] / distanceweight ) );
                }
              }
              else
              {
                if ( missing )
                {
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_ui[data_count++] = missingval_ui[dc];
                  }
                }
                else
                {
                  /* use the invalid positions array and set the data
                     value to zero */
                  if ( !DXSetElementInvalid( out_invalidhandle,
                                             invalid_count ) )
                    goto error;
                  for ( dc = 0; dc < shape[0]; dc++ )
                  {
                    new_data_ptr_ui[data_count++] = 0;
                  }
                }
              }
              break;
            default:
              DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
              goto error;
          }
          /* increment the invalid positions array pointer */
          invalid_count++;
        }
      }
    }

    if ( !DXChangedComponentValues( (Field)base, name ) )
      goto error;
    if ( !DXSetComponentValue( (Field)base, name, (Object)newcomponent ) )
      goto error;
    /* XXX and radius also (otherwise inv pos not nec) */
    /* should inv pos be done each time? (for each component) */
    if ( ( !missing ) && ( componentsdone == 1 ) )
    {
      if ( !DXSaveInvalidComponent( base, out_invalidhandle ) )
        goto error;
    }
    newcomponent = NULL;
    if ( !DXSetComponentAttribute( (Field)base, name, "dep",
                                   (Object)DXNewString( "positions" ) ) )
      goto error;

  component_done:
    compcount++;
    DXFree( (Pointer)missingval_f );
    missingval_f = NULL;
    DXFree( (Pointer)missingval_s );
    missingval_s = NULL;
    DXFree( (Pointer)missingval_i );
    missingval_i = NULL;
    DXFree( (Pointer)missingval_d );
    missingval_d = NULL;
    DXFree( (Pointer)missingval_b );
    missingval_b = NULL;
    DXFree( (Pointer)missingval_ui );
    missingval_ui = NULL;
    DXFree( (Pointer)missingval_ub );
    missingval_ub = NULL;
    DXFree( (Pointer)missingval_us );
    missingval_us = NULL;
    DXFree( (Pointer)bufdistance );
    bufdistance = NULL;
    DXFree( (Pointer)bufdata );
    bufdata = NULL;
    DXFree( (Pointer)data );
    data = NULL;
  }
  if ( componentsdone == 0 )
  {
    if ( missing )
    {
      DXWarning( "no position dependent components found; missing value given "
                 "cannot be used" );
    }
    else
    {
      /* need to allocate a buffer of length n to hold the nearest neighbors*/
      bufdistance = (float *)DXAllocateLocal( numnearest * sizeof( float ) );
      if ( !bufdistance )
        goto error;

      /* loop over all the grid points */
      inputposition_ptr = (float *)DXGetArrayData( positions_ino );
      for ( i = 0; i < counts[0]; i++ )
      {
        for ( j = 0; j < counts[1]; j++ )
        {
          for ( k = 0; k < counts[2]; k++ )
          {
            if ( regular )
            {
              switch ( shape_base[0] )
              {
                case( 3 ) :
                  gridposition = DXApply( DXVec( i, j, k ), matrix );
                  break;
                case( 2 ) :
                  gridposition2D = Apply2D( Vec2D( i, j ), matrix2D );
                  break;
                case( 1 ) :
                  gridposition1D = origin + i * delta;
                  break;
              }
            }
            else
            {
              switch ( shape_base[0] )
              {
                case( 3 ) :
                  dp = DXGetArrayEntry( handle, pos_count, scratch );
                  gridposition = DXPt( dp[0], dp[1], dp[2] );
                  break;
                case( 2 ) :
                  dp = DXGetArrayEntry( handle, pos_count, scratch );
                  gridposition2D.x = dp[0];
                  gridposition2D.y = dp[1];
                  break;
                case( 1 ) :
                  dp = DXGetArrayEntry( handle, pos_count, scratch );
                  gridposition1D = dp[0];
                  break;
              }
              pos_count++;
            }
            for ( ii = 0; ii < numnearest; ii++ )
            {
              bufdistance[ii] = DXD_MAX_FLOAT;
            }
            /* numvalid is how many valid data points are in the buffer
               for this particular grid point */
            numvalid = 0;
            for ( l = 0, tmpdata = old_data_ptr; l < numitemsino; l++ )
            {
              /* don't do anything if the input is invalid */
              if ( DXIsElementValid( in_invalidhandle, l ) )
              {
                switch ( shape_base[0] )
                {
                  case( 1 ) :
                    distance = ABS( gridposition1D - inputposition_ptr[l] );
                    break;
                  case( 2 ) :
                    distance = sqrt(
                        POW( gridposition2D.x - inputposition_ptr[2 * l], 2 ) +
                        POW( gridposition2D.y - inputposition_ptr[2 * l + 1],
                             2 ) );
                    break;
                  case( 3 ) :
                    distance = DXLength( DXSub(
                        gridposition, DXPt( inputposition_ptr[3 * l],
                                            inputposition_ptr[3 * l + 1],
                                            inputposition_ptr[3 * l + 2] ) ) );
                    break;
                }
                /* first try to insert into the already present buffer */
                for ( m = 0; m < numvalid; m++ )
                {
                  if ( distance < bufdistance[m] )
                  {
                    if ( numvalid < numnearest )
                      numvalid++;
                    for ( n = numvalid - 1; n > m; n-- )
                    {
                      bufdistance[n] = bufdistance[n - 1];
                    }
                    bufdistance[m] = distance;
                    goto filled_buf;
                  }
                }
                if ( numvalid == 0 )
                {
                  bufdistance[0] = distance;
                  numvalid++;
                }
              filled_buf:
                tmpdata += datastep;
                continue;
              }
            }

            /* have now looped over all the scattered points */

            if ( radius && ( bufdistance[0] > radiusvalue ) )
            {
              anyinvalid = 1;
              if ( !DXSetElementInvalid( out_invalidhandle, invalid_count ) )
                goto error;
            }
            /* increment the invalid positions array pointer */
            invalid_count++;
          }
        }
      }
      /* put the invalid component in the field */
      if ( anyinvalid )
      {
        if ( !DXSaveInvalidComponent( base, out_invalidhandle ) )
          goto error;
      }
      DXFree( (Pointer)bufdistance );
    }
  }

  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  DXFree( (Pointer)scratch );
  DXFreeArrayHandle( handle );
  if ( !DXEndField( base ) )
    goto error;
  return OK;

error:
  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  DXFree( (Pointer)scratch );
  DXFreeArrayHandle( handle );
  DXFree( (Pointer)bufdistance );
  DXFree( (Pointer)bufdata );
  DXFree( (Pointer)data );
  DXFree( (Pointer)missingval_f );
  DXFree( (Pointer)missingval_s );
  DXFree( (Pointer)missingval_i );
  DXFree( (Pointer)missingval_d );
  DXFree( (Pointer)missingval_b );
  DXFree( (Pointer)missingval_ui );
  DXFree( (Pointer)missingval_ub );
  DXFree( (Pointer)missingval_us );
  return ERROR;
}

static Error ConnectRadius( Field ino, Field base, float radius, float exponent,
                            Array missing )
{
  Array positions_base;

  if ( DXEmptyField( (Field)ino ) )
    return OK;

  positions_base = (Array)DXGetComponentValue( (Field)base, "positions" );
  if ( !positions_base )
  {
    DXSetError( ERROR_MISSING_DATA, "missing positions component in base" );
    goto error;
  }

  if ( DXQueryGridPositions( (Array)positions_base, NULL, NULL, NULL, NULL ) )
  {
    if ( !ConnectRadiusRegular( ino, base, radius, exponent, missing ) )
      goto error;
  }
  else
  {
    if ( !ConnectRadiusIrregular( ino, base, radius, exponent, missing ) )
      goto error;
  }
  if ( !DXEndField( base ) )
    goto error;
  return OK;
error:
  return ERROR;
}

static Error ConnectRadiusIrregular( Field ino, Field base, float radius,
                                     float exponent, Array missing )
{
  int componentsdone, valid;
  Array positions_base, positions_ino;
  Array oldcomponent = NULL, newcomponent = NULL;
  InvalidComponentHandle out_invalidhandle = NULL, in_invalidhandle = NULL;
  int numitemsbase, numitemsino;
  int rank_base, rank_ino, shape_base[30], shape_ino[30], i, j;
  int anyinvalid = 1, compcount, rank, shape[30];
  int numitems;
  Type type;
  ubyte *flagbuf = NULL;
  float *distancebuf = NULL;
  Category category;
  float *pos_ino_ptr, *pos_base_ptr;
  char *name;

  /* all points within radius of a given grid position will be averaged
   * (weighted) to give a data point to that position. If there are are
   * none within that range, the missing value is used. If missing is
   * NULL, then an invalid position indication is placed there and cull
   * will be used at the end */

  in_invalidhandle =
      DXCreateInvalidComponentHandle( (Object)ino, NULL, "positions" );
  if ( !in_invalidhandle )
    goto error;
  if ( !missing )
  {
    out_invalidhandle =
        DXCreateInvalidComponentHandle( (Object)base, NULL, "positions" );
    if ( !out_invalidhandle )
      goto error;
    anyinvalid = 0;
  }

  positions_base = (Array)DXGetComponentValue( (Field)base, "positions" );
  if ( !positions_base )
  {
    DXSetError( ERROR_MISSING_DATA, "base field has no positions" );
    goto error;
  }
  positions_ino = (Array)DXGetComponentValue( (Field)ino, "positions" );
  if ( !positions_ino )
  {
    DXSetError( ERROR_MISSING_DATA, "input field has no positions" );
    goto error;
  }
  DXGetArrayInfo( positions_base, &numitemsbase, &type, &category, &rank_base,
                  shape_base );
  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only real, floating point positions for base accepted" );
    goto error;
  }
  if ( ( rank_base > 1 ) || ( shape_base[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, or 3D positions for base accepted" );
    goto error;
  }
  DXGetArrayInfo( positions_ino, &numitemsino, &type, &category, &rank_ino,
                  shape_ino );
  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only real, floating point positions for input accepted" );
    goto error;
  }
  if ( ( rank_ino > 1 ) || ( shape_ino[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, or 3D positions for input accepted" );
    goto error;
  }

  if ( shape_base[0] != shape_ino[0] )
  {
    DXSetError( ERROR_DATA_INVALID,
                "base positions are %dD; input positions are %dD",
                shape_base[0], shape_ino[0] );
    goto error;
  }

  /* regular arrays? XXX */
  pos_ino_ptr = (float *)DXGetArrayData( positions_ino );

  /* make a buffer to hold the distances and flags */
  flagbuf = (ubyte *)DXAllocateLocalZero(
      ( numitemsbase * numitemsino * sizeof( ubyte ) ) );
  if ( !flagbuf )
    goto error;
  distancebuf = (float *)DXAllocateLocalZero(
      ( numitemsbase * numitemsino * sizeof( float ) ) );
  if ( !distancebuf )
    goto error;

  if ( out_invalidhandle )
  {
    DXSetAllValid( out_invalidhandle );
  }

  /* call irregular method */
  pos_base_ptr = (float *)DXGetArrayData( positions_base );
  if ( !( FillBuffersIrreg( numitemsbase, numitemsino, flagbuf, distancebuf,
                            pos_ino_ptr, pos_base_ptr, out_invalidhandle,
                            in_invalidhandle, &anyinvalid, shape_base[0],
                            radius ) ) )
    goto error;

  /* need to get the components we are going to copy into the new field */
  /* copied wholesale from gda */
  compcount = 0;
  componentsdone = 0;
  while ( NULL != ( oldcomponent = (Array)DXGetEnumeratedComponentValue(
                        (Field)ino, compcount, &name ) ) )
  {
    Object attr;

    if ( DXGetObjectClass( (Object)oldcomponent ) != CLASS_ARRAY )
      goto component_done;

    if ( !strcmp( name, "positions" ) || !strcmp( name, "connections" ) ||
         !strcmp( name, "invalid positions" ) )
      goto component_done;

    /* if the component refs anything, leave it */
    if ( DXGetComponentAttribute( (Field)ino, name, "ref" ) )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "der" );
    if ( attr )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "dep" );
    if ( !attr )
    {
      /* does it der? if so, throw it away */
      attr = DXGetComponentAttribute( (Field)ino, name, "der" );
      if ( attr )
        goto component_done;
      /* it doesn't ref, doesn't dep and doesn't der.
         Copy it to the output and quit */
      if ( !DXSetComponentValue( (Field)base, name, (Object)oldcomponent ) )
        goto error;
      goto component_done;
    }

    if ( DXGetObjectClass( attr ) != CLASS_STRING )
    {
      DXSetError( ERROR_DATA_INVALID, "dependency attribute not type string" );
      goto error;
    }

    if ( strcmp( DXGetString( (String)attr ), "positions" ) )
    {
      DXWarning( "Component '%s' is not dependent on positions!", name );
      DXWarning( "Regrid will remove '%s' and output the base grid 'data' "
                 "component if it exists",
                 name );
      goto component_done;
    }

    /* if the component happens to be "invalid positions" ignore it */
    if ( !strcmp( name, "invalid positions" ) )
      goto component_done;

    componentsdone++;

    DXGetArrayInfo( (Array)oldcomponent, &numitems, &type, &category, &rank,
                    shape );

    if ( rank == 0 )
      shape[0] = 1;
    if ( rank > 1 )
    {
      DXSetError( ERROR_NOT_IMPLEMENTED, "rank larger than 1 not implemented" );
      goto error;
    }
    if ( numitems != numitemsino )
    {
      DXSetError(
          ERROR_BAD_PARAMETER,
          "number of items in %s component not equal to number of positions",
          name );
      goto error;
    }
    newcomponent = DXNewArrayV( type, category, rank, shape );

    if ( !DXAddArrayData( newcomponent, 0, numitemsbase, NULL ) )
      goto error;

    /* here's where I use my preset buffers to fill the new array */
    /* XXX here do something about missing */
    if ( !FillDataArray( type, oldcomponent, newcomponent, flagbuf, distancebuf,
                         numitemsbase, numitemsino, shape[0], exponent, missing,
                         name ) )
      goto error;

    if ( !DXChangedComponentValues( (Field)base, name ) )
      goto error;
    if ( !DXSetComponentValue( (Field)base, name, (Object)newcomponent ) )
      goto error;
    newcomponent = NULL;
    if ( !DXSetComponentAttribute( (Field)base, name, "dep",
                                   (Object)DXNewString( "positions" ) ) )
      goto error;
  component_done:
    compcount++;
  }

  if ( componentsdone == 0 )
  {
    if ( missing )
    {
      DXWarning( "no position-dependent components found; missing value "
                 "provided cannot be used" );
    }
    else
    {
      for ( i = 0; i < numitemsbase; i++ )
      {
        valid = 0;
        for ( j = 0; j < numitemsino; j++ )
        {
          if ( flagbuf[i * numitemsino + j] )
          {
            /* there's at least one valid */
            valid = 1;
            break;
          }
        }
        if ( !valid )
        {
          DXSetElementInvalid( out_invalidhandle, i );
          anyinvalid = 1;
        }
      }
    }
  }
  else
  {
    if ( !missing )
    {
      for ( i = 0; i < numitemsbase; i++ )
      {
        valid = 0;
        for ( j = 0; j < numitemsino; j++ )
        {
          if ( flagbuf[i * numitemsino + j] )
          {
            /* there's at least one valid */
            valid = 1;
            break;
          }
        }
        if ( !valid )
        {
          DXSetElementInvalid( out_invalidhandle, i );
          anyinvalid = 1;
        }
      }
    }
  }

  if ( anyinvalid && out_invalidhandle )
  {
    if ( !DXSaveInvalidComponent( base, out_invalidhandle ) )
      goto error;
  }

  DXFree( (Pointer)flagbuf );
  DXFree( (Pointer)distancebuf );
  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  return OK;

error:
  DXDelete( (Object)newcomponent );
  DXFree( (Pointer)flagbuf );
  DXFree( (Pointer)distancebuf );
  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  return ERROR;
}

static Error ConnectRadiusRegular( Field ino, Field base, float radius,
                                   float exponent, Array missing )
{
  int componentsdone;
  Array positions_base, positions_ino;
  Array oldcomponent, newcomponent = NULL;
  InvalidComponentHandle in_invalidhandle = NULL, out_invalidhandle = NULL;
  int numitemsbase, numitemsino, cntr;
  int rank_base, rank_ino, shape_base[30], shape_ino[30], i, j, k, l, c;
  int compcount, rank, shape[30];
  float tmpval;
  int numitems, counts[8], knt;
  Type type;
  float xpos, ypos, zpos, xmin, ymin, zmin, xmax, ymax, zmax, pos;
  float min, max;
  Vector minvec, maxvec, posvec;
  Vector2D minvec2D, maxvec2D, posvec2D;
  Matrix matrix, invmatrix;
  Matrix2D matrix2D, invmatrix2D;
  Category category;
  float *pos_ino_ptr, distance, *weights_ptr = NULL;
  float *data_ptr = NULL;
  float rsquared;
  byte *exact_hit;
  Point box[8];
  char *name;
  float *data_in_f = NULL, *data_f = NULL, *missingval_f = NULL;
  byte *data_in_b = NULL, *data_b = NULL, *missingval_b = NULL;
  int *data_in_i = NULL, *data_i = NULL, *missingval_i = NULL;
  short *data_in_s = NULL, *data_s = NULL, *missingval_s = NULL;
  double *data_in_d = NULL, *data_d = NULL, *missingval_d = NULL;
  ubyte *data_in_ub = NULL, *data_ub = NULL, *missingval_ub = NULL;
  ushort *data_in_us = NULL, *data_us = NULL, *missingval_us = NULL;
  uint *data_in_ui = NULL, *data_ui = NULL, *missingval_ui = NULL;

  positions_base = (Array)DXGetComponentValue( (Field)base, "positions" );
  if ( !positions_base )
  {
    DXSetError( ERROR_MISSING_DATA, "base field has no positions" );
    goto error;
  }
  if ( !( DXQueryGridPositions( (Array)positions_base, NULL, counts,
                                (float *)&matrix.b, (float *)&matrix.A ) ) )
    goto error;
  /* all points within radius of a given grid position will be averaged
   * (weighted) to give a data point to that position. If there are are
   * none within that range, the missing value is used. If missing is
   * NULL, then an invalid position indication is placed there  */
  DXGetArrayInfo( positions_base, &numitemsbase, &type, &category, &rank_base,
                  shape_base );
  out_invalidhandle =
      DXCreateInvalidComponentHandle( (Object)base, NULL, "positions" );
  if ( !out_invalidhandle )
    goto error;

  DXSetAllInvalid( out_invalidhandle );

  /* if radius==-1, that's infinity. I really ought to handle it separately */
  /* XXX */
  if ( radius == -1 )
  {
    if ( !DXBoundingBox( (Object)ino, box ) )
      goto error;
    radius = 2 * DXLength( DXSub( box[0], box[7] ) );
  }

  positions_ino = (Array)DXGetComponentValue( (Field)ino, "positions" );
  if ( !positions_ino )
  {
    DXSetError( ERROR_MISSING_DATA, "input field has no positions" );
    goto error;
  }
  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only real, floating point positions for base accepted" );
    goto error;
  }
  if ( ( rank_base > 1 ) || ( shape_base[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, or 3D positions for base accepted" );
    goto error;
  }
  DXGetArrayInfo( positions_ino, &numitemsino, &type, &category, &rank_ino,
                  shape_ino );
  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only real, floating point positions for input accepted" );
    goto error;
  }
  if ( ( rank_ino > 1 ) || ( shape_ino[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, or 3D positions for input accepted" );
    goto error;
  }

  if ( shape_base[0] != shape_ino[0] )
  {
    DXSetError( ERROR_DATA_INVALID,
                "base positions are %dD; input positions are %dD",
                shape_base[0], shape_ino[0] );
    goto error;
  }

  /* regular arrays? XXX */
  pos_ino_ptr = (float *)DXGetArrayData( positions_ino );

  /* need to get the components we are going to copy into the new field */
  /* copied wholesale from gda */
  compcount = 0;
  componentsdone = 0;
  while ( NULL != ( oldcomponent = (Array)DXGetEnumeratedComponentValue(
                        (Field)ino, compcount, &name ) ) )
  {
    Object attr;

    if ( DXGetObjectClass( (Object)oldcomponent ) != CLASS_ARRAY )
      goto component_done;

    if ( !strcmp( name, "positions" ) || !strcmp( name, "connections" ) ||
         !strcmp( name, "invalid positions" ) )
      goto component_done;

    /* if the component refs anything, leave it */
    if ( DXGetComponentAttribute( (Field)ino, name, "ref" ) )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "der" );
    if ( attr )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "dep" );
    if ( !attr )
    {
      /* does it der? if so, throw it away */
      attr = DXGetComponentAttribute( (Field)ino, name, "der" );
      if ( attr )
        goto component_done;
      /* it doesn't ref, doesn't dep and doesn't der.
         Copy it to the output and quit */
      if ( !DXSetComponentValue( (Field)base, name, (Object)oldcomponent ) )
        goto error;
      goto component_done;
    }

    if ( DXGetObjectClass( attr ) != CLASS_STRING )
    {
      DXSetError( ERROR_DATA_INVALID, "dependency attribute not type string" );
      goto error;
    }

    if ( strcmp( DXGetString( (String)attr ), "positions" ) )
    {
      DXWarning( "Component '%s' is not dependent on positions!", name );
      DXWarning( "Regrid will remove '%s' and output the base grid 'data' "
                 "component if it exists",
                 name );
      goto component_done;
    }

    /* if the component happens to be "invalid positions" ignore it */
    if ( !strcmp( name, "invalid positions" ) )
      goto component_done;

    componentsdone++;

    DXGetArrayInfo( (Array)oldcomponent, &numitems, &type, &category, &rank,
                    shape );
    if ( rank == 0 )
      shape[0] = 1;

    /* check that the missing component, if present, matches the component
       we're working on */
    if ( missing )
      switch ( type )
      {
        case( TYPE_FLOAT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      float, type, shape[0], missingval_f );
          break;
        case( TYPE_INT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      int, type, shape[0], missingval_i );
          break;
        case( TYPE_BYTE ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      byte, type, shape[0], missingval_b );
          break;
        case( TYPE_DOUBLE ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      double, type, shape[0], missingval_d );
          break;
        case( TYPE_SHORT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      short, type, shape[0], missingval_s );
          break;
        case( TYPE_UBYTE ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      ubyte, type, shape[0], missingval_ub );
          break;
        case( TYPE_USHORT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      ushort, type, shape[0], missingval_us );
          break;
        case( TYPE_UINT ) :
          CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                      uint, type, shape[0], missingval_ui );
          break;
        default:
          DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
          goto error;
      }

    if ( rank == 0 )
      shape[0] = 1;
    if ( rank > 1 )
    {
      DXSetError( ERROR_NOT_IMPLEMENTED, "rank larger than 1 not implemented" );
      goto error;
    }
    if ( numitems != numitemsino )
    {
      DXSetError(
          ERROR_BAD_PARAMETER,
          "number of items in %s component not equal to number of positions",
          name );
      goto error;
    }

    weights_ptr =
        (float *)DXAllocateLocalZero( numitemsbase * sizeof( float ) );
    /* I need a floating point buffer to hold the weighted sum data */
    data_ptr = (float *)DXAllocateLocalZero( numitemsbase * sizeof( float ) *
                                             shape[0] );
    if ( !weights_ptr )
      goto error;
    if ( !data_ptr )
      goto error;

    /* new component is for the eventual output data */
    newcomponent = DXNewArrayV( type, category, rank, shape );
    if ( !DXAddArrayData( newcomponent, 0, numitemsbase, NULL ) )
      goto error;

    /* allocate the exact hit buf */
    exact_hit = (byte *)DXAllocateLocalZero( numitemsbase * sizeof( byte ) );

    switch ( type )
    {
      case( TYPE_FLOAT ) :
        data_in_f = (float *)DXGetArrayData( oldcomponent );
        data_f = (float *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_INT ) :
        data_in_i = (int *)DXGetArrayData( oldcomponent );
        data_i = (int *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_DOUBLE ) :
        data_in_d = (double *)DXGetArrayData( oldcomponent );
        data_d = (double *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_SHORT ) :
        data_in_s = (short *)DXGetArrayData( oldcomponent );
        data_s = (short *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_BYTE ) :
        data_in_b = (byte *)DXGetArrayData( oldcomponent );
        data_b = (byte *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_UBYTE ) :
        data_in_ub = (ubyte *)DXGetArrayData( oldcomponent );
        data_ub = (ubyte *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_USHORT ) :
        data_in_us = (ushort *)DXGetArrayData( oldcomponent );
        data_us = (ushort *)DXGetArrayData( newcomponent );
        break;
      case( TYPE_UINT ) :
        data_in_ui = (uint *)DXGetArrayData( oldcomponent );
        data_ui = (uint *)DXGetArrayData( newcomponent );
        break;
      default:
        DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
        goto error;
    }

#define INCREMENT_DATA_1D( type, data_in )                                     \
  {                                                                            \
    for ( i = 0; i < numitemsino; i++ )                                        \
    {                                                                          \
      xpos = pos_ino_ptr[i];                                                   \
      xmin = xpos - radius;                                                    \
      xmax = xpos + radius;                                                    \
      min = ( xmin / matrix.A[0][0] ) - ( matrix.b[0] / matrix.A[0][0] );      \
      max = ( xmax / matrix.A[0][0] ) - ( matrix.b[0] / matrix.A[0][0] );      \
      if ( min > max )                                                         \
      {                                                                        \
        tmpval = min;                                                          \
        min = max;                                                             \
        max = tmpval;                                                          \
      }                                                                        \
      /* minvec and maxvec are the indices corresponding to the                \
         bounding box around the sample point */                               \
      for ( j = CEIL( min ); j <= FLOOR( max ); j++ )                          \
        if ( j >= 0 && j < counts[0] )                                         \
        {                                                                      \
          /* here are the only ones which are in the bounding                  \
           * box */                                                            \
          pos = j * matrix.A[0][0] + matrix.b[0];                              \
          distance = ABS( pos - xpos );                                        \
          if ( distance <= radius )                                            \
          {                                                                    \
            DXSetElementValid( out_invalidhandle, j );                         \
            if ( ( distance > 0.0 ) && ( !exact_hit[j] ) )                     \
            {                                                                  \
              weights_ptr[j] += 1.0 / POW( distance, exponent );               \
              for ( c = 0; c < shape[0]; c++ )                                 \
                data_ptr[j * shape[0] + c] +=                                  \
                    (float)data_in[i * shape[0] + c] /                         \
                    POW( distance, exponent );                                 \
            }                                                                  \
            else if ( !exact_hit[j] )                                          \
            {                                                                  \
              /* we've hit right on; no more averaging here */                 \
              for ( c = 0; c < shape[0]; c++ )                                 \
                data_ptr[j * shape[0] + c] = (float)data_in[i * shape[0] + c]; \
              weights_ptr[j] = 1;                                              \
              exact_hit[j] = 1;                                                \
            }                                                                  \
          }                                                                    \
        }                                                                      \
    }                                                                          \
  }

#define INCREMENT_DATA_2D( type, data_in )                                   \
  {                                                                          \
    Matrix2D tmpmat;                                                         \
    matrix2D.A[0][0] = matrix.A[0][0];                                       \
    matrix2D.A[1][0] = matrix.A[0][1];                                       \
    matrix2D.A[0][1] = matrix.A[0][2];                                       \
    matrix2D.A[1][1] = matrix.A[1][0];                                       \
    matrix2D.b[0] = matrix.b[0];                                             \
    matrix2D.b[1] = matrix.b[1];                                             \
    invmatrix2D = Invert2D( matrix2D );                                      \
    /* loop over all scattered points */                                     \
    rsquared = radius * radius;                                              \
    for ( i = 0; i < numitemsino; i++ )                                      \
    {                                                                        \
      xpos = pos_ino_ptr[2 * i];                                             \
      ypos = pos_ino_ptr[2 * i + 1];                                         \
      xmin = xpos - radius;                                                  \
      xmax = xpos + radius;                                                  \
      ymin = ypos - radius;                                                  \
      ymax = ypos + radius;                                                  \
      tmpmat = invmatrix2D;                                                  \
      minvec2D = Apply2D( Vec2D( xmin, ymin ), tmpmat );                     \
      maxvec2D = Apply2D( Vec2D( xmax, ymax ), tmpmat );                     \
      if ( minvec2D.x > maxvec2D.x )                                         \
      {                                                                      \
        tmpval = minvec2D.x;                                                 \
        minvec2D.x = maxvec2D.x;                                             \
        maxvec2D.x = tmpval;                                                 \
      }                                                                      \
      if ( minvec2D.y > maxvec2D.y )                                         \
      {                                                                      \
        tmpval = minvec2D.y;                                                 \
        minvec2D.y = maxvec2D.y;                                             \
        maxvec2D.y = tmpval;                                                 \
      }                                                                      \
      /* minvec and maxvec are the indices corresponding to                  \
         the bounding box around the sample point */                         \
      tmpmat = matrix2D;                                                     \
      for ( j = CEIL( minvec2D.x ); j <= FLOOR( maxvec2D.x ); j++ )          \
        if ( j >= 0 && j < counts[0] )                                       \
          for ( k = CEIL( minvec2D.y ); k <= FLOOR( maxvec2D.y ); k++ )      \
            if ( k >= 0 && k < counts[1] )                                   \
            {                                                                \
              /* here are the only ones which are in the                     \
                 bounding box */                                             \
              posvec2D = Apply2D( Vec2D( j, k ), tmpmat );                   \
              distance = ( posvec2D.x - xpos ) * ( posvec2D.x - xpos ) +     \
                         ( posvec2D.y - ypos ) * ( posvec2D.y - ypos );      \
              if ( distance <= rsquared )                                    \
              {                                                              \
                distance = sqrt( distance );                                 \
                cntr = j * counts[1] + k;                                    \
                DXSetElementValid( out_invalidhandle, cntr );                \
                if ( ( distance > 0.0 ) && ( !exact_hit[cntr] ) )            \
                {                                                            \
                  weights_ptr[cntr] += 1.0 / POW( distance, exponent );      \
                  for ( c = 0; c < shape[0]; c++ )                           \
                    data_ptr[j * counts[1] * shape[0] + k * shape[0] + c] += \
                        (float)data_in[i * shape[0] + c] /                   \
                        POW( distance, exponent );                           \
                }                                                            \
                else if ( !exact_hit[cntr] )                                 \
                {                                                            \
                  for ( c = 0; c < shape[0]; c++ )                           \
                    data_ptr[j * counts[1] * shape[0] + k * shape[0]] =      \
                        (float)data_in[i * shape[0] + c];                    \
                  weights_ptr[cntr] = 1;                                     \
                  exact_hit[cntr] = 1;                                       \
                }                                                            \
              }                                                              \
            }                                                                \
    }                                                                        \
  }

#define INCREMENT_DATA_3D( type, data_in )                                   \
  {                                                                          \
    Matrix tmpmat;                                                           \
    invmatrix = DXInvert( matrix );                                          \
    /* loop over all scattered points */                                     \
    for ( knt = 0, i = 0; i < numitemsino; i++ )                             \
    {                                                                        \
      xpos = pos_ino_ptr[knt];                                               \
      ypos = pos_ino_ptr[knt + 1];                                           \
      zpos = pos_ino_ptr[knt + 2];                                           \
      knt += 3;                                                              \
      xmin = xpos - radius;                                                  \
      xmax = xpos + radius;                                                  \
      ymin = ypos - radius;                                                  \
      ymax = ypos + radius;                                                  \
      zmin = zpos - radius;                                                  \
      zmax = zpos + radius;                                                  \
      tmpmat = invmatrix;                                                    \
      minvec = DXApply( DXVec( xmin, ymin, zmin ), tmpmat );                 \
      maxvec = DXApply( DXVec( xmax, ymax, zmax ), tmpmat );                 \
      if ( minvec.x > maxvec.x )                                             \
      {                                                                      \
        tmpval = minvec.x;                                                   \
        minvec.x = maxvec.x;                                                 \
        maxvec.x = tmpval;                                                   \
      }                                                                      \
      if ( minvec.y > maxvec.y )                                             \
      {                                                                      \
        tmpval = minvec.y;                                                   \
        minvec.y = maxvec.y;                                                 \
        maxvec.y = tmpval;                                                   \
      }                                                                      \
      if ( minvec.z > maxvec.z )                                             \
      {                                                                      \
        tmpval = minvec.z;                                                   \
        minvec.z = maxvec.z;                                                 \
        maxvec.z = tmpval;                                                   \
      }                                                                      \
      /* minvec and maxvec are the indices corresponding to the */           \
      /* bounding box around the sample point */                             \
      tmpmat = matrix;                                                       \
      rsquared = radius * radius;                                            \
      for ( j = CEIL( minvec.x ); j <= FLOOR( maxvec.x ); j++ )              \
        if ( j >= 0 && j < counts[0] )                                       \
          for ( k = CEIL( minvec.y ); k <= FLOOR( maxvec.y ); k++ )          \
            if ( k >= 0 && k < counts[1] )                                   \
              for ( l = CEIL( minvec.z ); l <= FLOOR( maxvec.z ); l++ )      \
                if ( l >= 0 && l < counts[2] )                               \
                {                                                            \
                  /* here are the only ones which are in the                 \
                     bounding box */                                         \
                  posvec = DXApply( DXVec( j, k, l ), tmpmat );              \
                  distance = ( posvec.x - xpos ) * ( posvec.x - xpos ) +     \
                             ( posvec.y - ypos ) * ( posvec.y - ypos ) +     \
                             ( posvec.z - zpos ) * ( posvec.z - zpos );      \
                  if ( distance <= rsquared )                                \
                  {                                                          \
                    distance = sqrt( distance );                             \
                    cntr = j * counts[2] * counts[1] + k * counts[2] + l;    \
                    DXSetElementValid( out_invalidhandle, cntr );            \
                    if ( ( distance > 0.0 ) && ( !exact_hit[cntr] ) )        \
                    {                                                        \
                      weights_ptr[cntr] += 1.0 / POW( distance, exponent );  \
                      for ( c = 0; c < shape[0]; c++ )                       \
                        data_ptr                                             \
                            [j * counts[2] * counts[1] * shape[0] +          \
                             k * counts[2] * shape[0] + l * shape[0] + c] += \
                            (float)data_in[i * shape[0] + c] /               \
                            POW( distance, exponent );                       \
                    }                                                        \
                    else if ( !exact_hit[cntr] )                             \
                    {                                                        \
                      for ( c = 0; c < shape[0]; c++ )                       \
                        data_ptr                                             \
                            [j * counts[2] * counts[1] * shape[0] +          \
                             k * counts[2] * shape[0] + l * shape[0] + c] =  \
                                (float)data_in[i * shape[0] + c];            \
                      weights_ptr[cntr] = 1;                                 \
                      exact_hit[cntr] = 1;                                   \
                    }                                                        \
                  }                                                          \
                }                                                            \
    }                                                                        \
  }

    if ( shape_base[0] == 1 )
    {
      switch ( type )
      {
        case( TYPE_FLOAT ) :
          INCREMENT_DATA_1D( float, data_in_f );
          break;
        case( TYPE_INT ) :
          INCREMENT_DATA_1D( int, data_in_i );
          break;
        case( TYPE_BYTE ) :
          INCREMENT_DATA_1D( byte, data_in_b );
          break;
        case( TYPE_SHORT ) :
          INCREMENT_DATA_1D( short, data_in_s );
          break;
        case( TYPE_DOUBLE ) :
          INCREMENT_DATA_1D( double, data_in_d );
          break;
        case( TYPE_UBYTE ) :
          INCREMENT_DATA_1D( ubyte, data_in_ub );
          break;
        case( TYPE_UINT ) :
          INCREMENT_DATA_1D( uint, data_in_ui );
          break;
        case( TYPE_USHORT ) :
          INCREMENT_DATA_1D( ushort, data_in_us );
          break;
        default:
          DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
          goto error;
      }
    }
    else if ( shape_base[0] == 2 )
    {
      switch ( type )
      {
        case( TYPE_FLOAT ) :
          INCREMENT_DATA_2D( float, data_in_f );
          break;
        case( TYPE_INT ) :
          INCREMENT_DATA_2D( int, data_in_i );
          break;
        case( TYPE_BYTE ) :
          INCREMENT_DATA_2D( byte, data_in_b );
          break;
        case( TYPE_SHORT ) :
          INCREMENT_DATA_2D( short, data_in_s );
          break;
        case( TYPE_DOUBLE ) :
          INCREMENT_DATA_2D( double, data_in_d );
          break;
        case( TYPE_UBYTE ) :
          INCREMENT_DATA_2D( ubyte, data_in_ub );
          break;
        case( TYPE_UINT ) :
          INCREMENT_DATA_2D( uint, data_in_ui );
          break;
        case( TYPE_USHORT ) :
          INCREMENT_DATA_2D( ushort, data_in_us );
          break;
        default:
          DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
          goto error;
      }
    }
    else if ( shape_base[0] == 3 )
    {
      switch ( type )
      {
        case( TYPE_FLOAT ) :
          INCREMENT_DATA_3D( float, data_in_f );
          break;
        case( TYPE_INT ) :
          INCREMENT_DATA_3D( int, data_in_i );
          break;
        case( TYPE_BYTE ) :
          INCREMENT_DATA_3D( byte, data_in_b );
          break;
        case( TYPE_SHORT ) :
          INCREMENT_DATA_3D( short, data_in_s );
          break;
        case( TYPE_DOUBLE ) :
          INCREMENT_DATA_3D( double, data_in_d );
          break;
        case( TYPE_UBYTE ) :
          INCREMENT_DATA_3D( ubyte, data_in_ub );
          break;
        case( TYPE_UINT ) :
          INCREMENT_DATA_3D( uint, data_in_ui );
          break;
        case( TYPE_USHORT ) :
          INCREMENT_DATA_3D( ushort, data_in_us );
          break;
        default:
          DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
          goto error;
      }
    }
    else
    {
      DXSetError( ERROR_DATA_INVALID,
                  "only 1D, 2D, and 3D input positions accepted" );
      goto error;
    }

    /* now fill the actual data array (int or float or whatever)
       with the computed sums */
    switch ( type )
    {
      case( TYPE_FLOAT ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_f[i * shape[0] + c] =
                  data_ptr[i * shape[0] + c] / weights_ptr[i];
          else
          {
            /* XXX here do something about missing */
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_f[i * shape[0] + c] = missingval_f[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_f[i * shape[0] + c] = 0;
          }
        }
        break;
      case( TYPE_INT ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_i[i * shape[0] + c] =
                  (int)( RND( data_ptr[i * shape[0] + c] / weights_ptr[i] ) );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_i[i * shape[0] + c] = (int)missingval_i[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_i[i * shape[0] + c] = (int)0;
          }
        }
        break;
      case( TYPE_SHORT ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_s[i * shape[0] + c] =
                  (short)( RND( data_ptr[i * shape[0] + c] / weights_ptr[i] ) );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_s[i * shape[0] + c] = (short)missingval_s[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_s[i * shape[0] + c] = (short)0;
          }
        }
        break;
      case( TYPE_BYTE ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_b[i * shape[0] + c] = ( byte )(
                  RND( data_ptr[i * shape[0] + c] / weights_ptr[i] ) );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_b[i * shape[0] + c] = (byte)missingval_b[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_b[i * shape[0] + c] = (byte)0;
          }
        }
        break;
      case( TYPE_DOUBLE ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_d[i * shape[0] + c] =
                  (double)( data_ptr[i * shape[0] + c] / weights_ptr[i] );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_d[i * shape[0] + c] = (double)missingval_d[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_d[i * shape[0] + c] = (double)0;
          }
        }
        break;
      case( TYPE_UINT ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_ui[i * shape[0] + c] = ( uint )(
                  RND( data_ptr[i * shape[0] + c] / weights_ptr[i] ) );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_ui[i * shape[0] + c] = (uint)missingval_ui[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_ui[i * shape[0] + c] = (uint)0;
          }
        }
        break;
      case( TYPE_USHORT ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_us[i * shape[0] + c] = ( ushort )(
                  RND( data_ptr[i * shape[0] + c] / weights_ptr[i] ) );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_us[i * shape[0] + c] = (ushort)missingval_us[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_us[i * shape[0] + c] = (ushort)0;
          }
        }
        break;
      case( TYPE_UBYTE ) :
        for ( i = 0; i < numitemsbase; i++ )
        {
          if ( DXIsElementValid( out_invalidhandle, i ) )
            for ( c = 0; c < shape[0]; c++ )
              data_ub[i * shape[0] + c] = ( ubyte )(
                  RND( data_ptr[i * shape[0] + c] / weights_ptr[i] ) );
          else
          {
            if ( missing )
              for ( c = 0; c < shape[0]; c++ )
                data_ub[i * shape[0] + c] = (ubyte)missingval_ub[c];
            else
              for ( c = 0; c < shape[0]; c++ )
                data_ub[i * shape[0] + c] = (ubyte)0;
          }
        }
        break;
      default:
        DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
        goto error;
    }

    if ( !DXSetComponentValue( (Field)base, name, (Object)newcomponent ) )
      goto error;
    if ( !DXSetComponentAttribute( (Field)base, name, "dep",
                                   (Object)DXNewString( "positions" ) ) )
      goto error;
    newcomponent = NULL;
    if ( !DXChangedComponentValues( (Field)base, name ) )
      goto error;
  component_done:
    compcount++;
    DXFree( (Pointer)missingval_f );
    missingval_f = NULL;
    DXFree( (Pointer)missingval_b );
    missingval_b = NULL;
    DXFree( (Pointer)missingval_i );
    missingval_i = NULL;
    DXFree( (Pointer)missingval_s );
    missingval_s = NULL;
    DXFree( (Pointer)missingval_d );
    missingval_d = NULL;
    DXFree( (Pointer)missingval_ub );
    missingval_ub = NULL;
    DXFree( (Pointer)missingval_us );
    missingval_us = NULL;
    DXFree( (Pointer)missingval_ui );
    missingval_ui = NULL;
  }

  if ( componentsdone == 0 )
  {
    if ( missing )
    {
      DXWarning( "no position-dependent components found; cannot use missing "
                 "value provided" );
    }
    else
    {
      switch ( shape_base[0] )
      {
        case( 1 ) :
          for ( i = 0; i < numitemsino; i++ )
          {
            xpos = pos_ino_ptr[i];
            xmin = xpos - radius;
            xmax = xpos + radius;
            min = ( xmin / matrix.A[0][0] ) - ( matrix.b[0] / matrix.A[0][0] );
            max = ( xmax / matrix.A[0][0] ) - ( matrix.b[0] / matrix.A[0][0] );
            /* minvec and maxvec are the indices corresponding to the
               bounding box around the sample point */
            for ( j = CEIL( min ); j <= FLOOR( max ); j++ )
              if ( j >= 0 && j < counts[0] )
              {
                /* here are the only ones which are in the bounding
                 * box */
                pos = j * matrix.A[0][0] + matrix.b[0];
                distance = pos - xpos;
                if ( distance <= radius )
                {
                  DXSetElementValid( out_invalidhandle, j );
                }
              }
          }
          break;
        case( 2 ) :
        {
          Matrix2D tmpmat;
          matrix2D.A[0][0] = matrix.A[0][0];
          matrix2D.A[1][0] = matrix.A[0][1];
          matrix2D.A[0][1] = matrix.A[0][2];
          matrix2D.A[1][1] = matrix.A[1][0];
          matrix2D.b[0] = matrix.b[0];
          matrix2D.b[1] = matrix.b[1];
          invmatrix2D = Invert2D( matrix2D );
          /* loop over all scattered points */
          for ( i = 0; i < numitemsino; i++ )
          {
            xpos = pos_ino_ptr[2 * i];
            ypos = pos_ino_ptr[2 * i + 1];
            xmin = xpos - radius;
            xmax = xpos + radius;
            ymin = ypos - radius;
            ymax = ypos + radius;
            tmpmat = invmatrix2D;
            minvec2D = Apply2D( Vec2D( xmin, ymin ), tmpmat );
            maxvec2D = Apply2D( Vec2D( xmax, ymax ), tmpmat );
            /* minvec and maxvec are the indices corresponding to
               the bounding box around the sample point */
            tmpmat = matrix2D;
            for ( j = CEIL( minvec2D.x ); j <= FLOOR( maxvec2D.x ); j++ )
              if ( j >= 0 && j < counts[0] )
                for ( k = CEIL( minvec2D.y ); k <= FLOOR( maxvec2D.y ); k++ )
                  if ( k >= 0 && k < counts[1] )
                  {
                    /* here are the only ones which are in the
                       bounding box */
                    posvec2D = Apply2D( Vec2D( j, k ), tmpmat );
                    distance = POW( posvec2D.x - xpos, 2 ) +
                               POW( posvec2D.y - ypos, 2 );
                    distance = sqrt( distance );
                    if ( distance <= radius )
                    {
                      cntr = j * counts[1] + k;
                      DXSetElementValid( out_invalidhandle, cntr );
                    }
                  }
          }
          break;
        }
        case( 3 ) :
        {
          Matrix tmpmat;
          invmatrix = DXInvert( matrix );
          /* loop over all scattered points */
          for ( i = 0; i < numitemsino; i++ )
          {
            xpos = pos_ino_ptr[3 * i];
            ypos = pos_ino_ptr[3 * i + 1];
            zpos = pos_ino_ptr[3 * i + 2];
            xmin = xpos - radius;
            xmax = xpos + radius;
            ymin = ypos - radius;
            ymax = ypos + radius;
            zmin = zpos - radius;
            zmax = zpos + radius;
            tmpmat = invmatrix;
            minvec = DXApply( DXVec( xmin, ymin, zmin ), tmpmat );
            maxvec = DXApply( DXVec( xmax, ymax, zmax ), tmpmat );
            /* minvec and maxvec are the indices corresponding to the */
            /* bounding box around the sample point */
            tmpmat = matrix;
            for ( j = CEIL( minvec.x ); j <= FLOOR( maxvec.x ); j++ )
              if ( j >= 0 && j < counts[0] )
                for ( k = CEIL( minvec.y ); k <= FLOOR( maxvec.y ); k++ )
                  if ( k >= 0 && k < counts[1] )
                    for ( l = CEIL( minvec.z ); l <= FLOOR( maxvec.z ); l++ )
                      if ( l >= 0 && l < counts[2] )
                      {
                        /* here are the only ones which are in the
                           bounding box */
                        posvec = DXApply( DXVec( j, k, l ), tmpmat );
                        distance = POW( posvec.x - xpos, 2 ) +
                                   POW( posvec.y - ypos, 2 ) +
                                   POW( posvec.z - zpos, 2 );
                        distance = sqrt( distance );
                        if ( distance <= radius )
                        {
                          cntr = j * counts[2] * counts[1] + k * counts[2] + l;
                          DXSetElementValid( out_invalidhandle, cntr );
                        }
                      }
          }
          break;
        }
        default:
          DXSetError( ERROR_DATA_INVALID,
                      "grid must be 1-, 2-, or 3-dimensional" );
          goto error;
          break;
      }
    }
  }

  if ( !missing )
  {
    if ( !DXSaveInvalidComponent( base, out_invalidhandle ) )
      goto error;
  }

  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  DXFree( (Pointer)weights_ptr );
  return OK;

error:
  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  DXDelete( (Object)newcomponent );
  DXFree( (Pointer)weights_ptr );
  DXFree( (Pointer)missingval_f );
  DXFree( (Pointer)missingval_b );
  DXFree( (Pointer)missingval_i );
  DXFree( (Pointer)missingval_s );
  DXFree( (Pointer)missingval_d );
  DXFree( (Pointer)missingval_ub );
  DXFree( (Pointer)missingval_us );
  DXFree( (Pointer)missingval_ui );
  return ERROR;
}

static Error FillDataArray( Type type, Array oldcomponent, Array newcomponent,
                            ubyte *flagbuf, float *distancebuf,
                            int numitemsbase, int numitemsino, int shape,
                            float exponent, Array missing, char *name )
{
  float weight, distance;
  int i, j, k;
  float *buffer = NULL;
  float *d_in_f = NULL, *d_f = NULL, *missingval_f = NULL;
  int *d_in_i = NULL, *d_i = NULL, *missingval_i = NULL;
  short *d_in_s = NULL, *d_s = NULL, *missingval_s = NULL;
  double *d_in_d = NULL, *d_d = NULL, *missingval_d = NULL;
  byte *d_in_b = NULL, *d_b = NULL, *missingval_b = NULL;
  uint *d_in_ui = NULL, *d_ui = NULL, *missingval_ui = NULL;
  ushort *d_in_us = NULL, *d_us = NULL, *missingval_us = NULL;
  ubyte *d_in_ub = NULL, *d_ub = NULL, *missingval_ub = NULL;

  buffer = (float *)DXAllocateLocal( shape * sizeof( float ) );
  if ( !buffer )
    goto error;

  /* check that the missing component, if present, matches the component
     we're working on */
  if ( missing )
    switch ( type )
    {
      case( TYPE_FLOAT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    float, type, shape, missingval_f );
        break;
      case( TYPE_INT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ), int,
                                    type, shape, missingval_i );
        break;
      case( TYPE_BYTE ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    byte, type, shape, missingval_b );
        break;
      case( TYPE_DOUBLE ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    double, type, shape, missingval_d );
        break;
      case( TYPE_SHORT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    short, type, shape, missingval_s );
        break;
      case( TYPE_UBYTE ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    ubyte, type, shape, missingval_ub );
        break;
      case( TYPE_USHORT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    ushort, type, shape, missingval_us );
        break;
      case( TYPE_UINT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    uint, type, shape, missingval_ui );
        break;
      default:
        DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
        goto error;
    }

  switch ( type )
  {
    case( TYPE_FLOAT ) :
      d_in_f = (float *)DXGetArrayData( oldcomponent );
      d_f = (float *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_INT ) :
      d_in_i = (int *)DXGetArrayData( oldcomponent );
      d_i = (int *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_SHORT ) :
      d_in_s = (short *)DXGetArrayData( oldcomponent );
      d_s = (short *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_DOUBLE ) :
      d_in_d = (double *)DXGetArrayData( oldcomponent );
      d_d = (double *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_BYTE ) :
      d_in_b = (byte *)DXGetArrayData( oldcomponent );
      d_b = (byte *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_UBYTE ) :
      d_in_ub = (ubyte *)DXGetArrayData( oldcomponent );
      d_ub = (ubyte *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_USHORT ) :
      d_in_us = (ushort *)DXGetArrayData( oldcomponent );
      d_us = (ushort *)DXGetArrayData( newcomponent );
      break;
    case( TYPE_UINT ) :
      d_in_ui = (uint *)DXGetArrayData( oldcomponent );
      d_ui = (uint *)DXGetArrayData( newcomponent );
      break;
    default:
      DXSetError( ERROR_DATA_INVALID, "unrecognized data type" );
      goto error;
  }

  for ( i = 0; i < numitemsbase; i++ )
  {
    weight = 0;
    for ( k = 0; k < shape; k++ )
    {
      buffer[k] = 0.0;
    }
    for ( j = 0; j < numitemsino; j++ )
    {
      if ( flagbuf[i * numitemsino + j] )
      {
        distance = distancebuf[i * numitemsino + j];
        if ( distance != 0 )
        {
          for ( k = 0; k < shape; k++ )
          {
            switch ( type )
            {
              case( TYPE_FLOAT ) :
                buffer[k] +=
                    (float)d_in_f[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_INT ) :
                buffer[k] +=
                    (float)d_in_i[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_SHORT ) :
                buffer[k] +=
                    (float)d_in_s[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_BYTE ) :
                buffer[k] +=
                    (float)d_in_b[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_DOUBLE ) :
                buffer[k] +=
                    (float)d_in_d[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_UBYTE ) :
                buffer[k] +=
                    (float)d_in_ub[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_USHORT ) :
                buffer[k] +=
                    (float)d_in_us[j * shape + k] / POW( distance, exponent );
                break;
              case( TYPE_UINT ) :
                buffer[k] +=
                    (float)d_in_ui[j * shape + k] / POW( distance, exponent );
                break;
              default:
                DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
                goto error;
            }
          }
          weight += 1 / POW( distance, exponent );
        }
        else
        {
          /* if distance == 0, then simply use the data value */
          for ( k = 0; k < shape; k++ )
          {
            switch ( type )
            {
              case( TYPE_FLOAT ) :
                buffer[k] = (float)d_in_f[j];
                break;
              case( TYPE_INT ) :
                buffer[k] = (float)d_in_i[j];
                break;
              case( TYPE_BYTE ) :
                buffer[k] = (float)d_in_b[j];
                break;
              case( TYPE_DOUBLE ) :
                buffer[k] = (float)d_in_d[j];
                break;
              case( TYPE_SHORT ) :
                buffer[k] = (float)d_in_s[j];
                break;
              case( TYPE_UINT ) :
                buffer[k] = (float)d_in_ui[j];
                break;
              case( TYPE_USHORT ) :
                buffer[k] = (float)d_in_us[j];
                break;
              case( TYPE_UBYTE ) :
                buffer[k] = (float)d_in_us[j];
                break;
              default:
                DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
                goto error;
            }
          }
          weight = 1;
          continue;
        }
      }
    }
    if ( weight > 0 )
    {
      for ( k = 0; k < shape; k++ )
      {
        switch ( type )
        {
          case( TYPE_FLOAT ) :
            d_f[shape * i + k] = (float)( buffer[k] / weight );
            break;
          case( TYPE_DOUBLE ) :
            d_d[shape * i + k] = (double)( buffer[k] / weight );
            break;
          case( TYPE_INT ) :
            d_i[shape * i + k] = (int)( RND( buffer[k] / weight ) );
            break;
          case( TYPE_SHORT ) :
            d_s[shape * i + k] = (short)( RND( buffer[k] / weight ) );
            break;
          case( TYPE_BYTE ) :
            d_b[shape * i + k] = ( byte )( RND( buffer[k] / weight ) );
            break;
          case( TYPE_UBYTE ) :
            d_ub[shape * i + k] = ( ubyte )( RND( buffer[k] / weight ) );
            break;
          case( TYPE_UINT ) :
            d_ui[shape * i + k] = ( uint )( RND( buffer[k] / weight ) );
            break;
          case( TYPE_USHORT ) :
            d_us[shape * i + k] = ( ushort )( RND( buffer[k] / weight ) );
            break;
          default:
            DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
            goto error;
        }
      }
    }
    else
    {
      if ( missing )
        for ( k = 0; k < shape; k++ )
        {
          switch ( type )
          {
            case( TYPE_FLOAT ) :
              d_f[i * shape + k] = (float)missingval_f[k];
              break;
            case( TYPE_INT ) :
              d_i[i * shape + k] = (int)missingval_i[k];
              break;
            case( TYPE_BYTE ) :
              d_b[i * shape + k] = (byte)missingval_b[k];
              break;
            case( TYPE_DOUBLE ) :
              d_d[i * shape + k] = (double)missingval_d[k];
              break;
            case( TYPE_SHORT ) :
              d_s[i * shape + k] = (short)missingval_s[k];
              break;
            case( TYPE_UBYTE ) :
              d_ub[i * shape + k] = (ubyte)missingval_ub[k];
              break;
            case( TYPE_UINT ) :
              d_ui[i * shape + k] = (uint)missingval_ui[k];
              break;
            case( TYPE_USHORT ) :
              d_us[i * shape + k] = (ushort)missingval_us[k];
              break;
            default:
              DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
              goto error;
          }
        }
      else
        for ( k = 0; k < shape; k++ )
        {
          switch ( type )
          {
            case( TYPE_FLOAT ) :
              d_f[i * shape + k] = 0.0;
              break;
            case( TYPE_INT ) :
              d_i[i * shape + k] = 0.0;
              break;
            case( TYPE_BYTE ) :
              d_b[i * shape + k] = 0.0;
              break;
            case( TYPE_DOUBLE ) :
              d_d[i * shape + k] = 0.0;
              break;
            case( TYPE_SHORT ) :
              d_s[i * shape + k] = 0.0;
              break;
            case( TYPE_UBYTE ) :
              d_ub[i * shape + k] = 0.0;
              break;
            case( TYPE_UINT ) :
              d_ui[i * shape + k] = 0.0;
              break;
            case( TYPE_USHORT ) :
              d_us[i * shape + k] = 0.0;
              break;
            default:
              DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
              goto error;
          }
        }
    }
  }
  DXFree( (Pointer)missingval_f );
  missingval_f = NULL;
  DXFree( (Pointer)missingval_i );
  missingval_i = NULL;
  DXFree( (Pointer)missingval_d );
  missingval_d = NULL;
  DXFree( (Pointer)missingval_s );
  missingval_s = NULL;
  DXFree( (Pointer)missingval_b );
  missingval_b = NULL;
  DXFree( (Pointer)missingval_ub );
  missingval_ub = NULL;
  DXFree( (Pointer)missingval_ui );
  missingval_ui = NULL;
  DXFree( (Pointer)missingval_us );
  missingval_us = NULL;
  return OK;
error:
  DXFree( (Pointer)buffer );
  DXFree( (Pointer)missingval_f );
  DXFree( (Pointer)missingval_i );
  DXFree( (Pointer)missingval_d );
  DXFree( (Pointer)missingval_s );
  DXFree( (Pointer)missingval_b );
  DXFree( (Pointer)missingval_ub );
  DXFree( (Pointer)missingval_ui );
  DXFree( (Pointer)missingval_us );
  return ERROR;
}

static Error FillBuffersIrreg( int numitemsbase, int numitemsino,
                               ubyte *flagbuf, float *distancebuf,
                               float *pos_ino_ptr, float *pos_base_ptr,
                               InvalidComponentHandle out_invalidhandle,
                               InvalidComponentHandle in_invalidhandle,
                               int *anyinvalid, int shape_base, float radius )
{
  int i, j, anyvalid;
  float distance = 0;

  /* there's a special check in here for radius == -1 (infinity) */

  /* fill the buffers */
  if ( shape_base == 1 )
  {
    for ( i = 0; i < numitemsbase; i++ )
    {
      anyvalid = 0;
      for ( j = 0; j < numitemsino; j++ )
      {
        if ( DXIsElementValid( in_invalidhandle, j ) )
        {
          if ( radius != -1 )
          {
            distance = ABS( pos_ino_ptr[j] - pos_base_ptr[i] );
            if ( distance <= radius )
            {
              anyvalid = 1;
              flagbuf[i * numitemsino + j] = 1;
              distancebuf[i * numitemsino + j] = distance;
            }
          }
          else
          {
            anyvalid = 1;
            flagbuf[i * numitemsino + j] = 1;
            distancebuf[i * numitemsino + j] = distance;
          }
        }
      }
      if ( !anyvalid )
      {
        if ( out_invalidhandle )
          DXSetElementInvalid( out_invalidhandle, i );
        *anyinvalid = 1;
      }
    }
  }

  else if ( shape_base == 2 )
  {
    for ( i = 0; i < numitemsbase; i++ )
    {
      anyvalid = 0;
      for ( j = 0; j < numitemsino; j++ )
      {
        if ( DXIsElementValid( in_invalidhandle, j ) )
        {
          if ( radius != -1 )
          {
            if ( ( ( pos_ino_ptr[2 * j] - pos_base_ptr[2 * i] ) < radius ) &&
                 ( ( pos_ino_ptr[2 * j + 1] - pos_base_ptr[2 * i + 1] ) <
                   radius ) )
            {
              distance =
                  POW( pos_ino_ptr[2 * j] - pos_base_ptr[2 * i], 2 ) +
                  POW( pos_ino_ptr[2 * j + 1] - pos_base_ptr[2 * i + 1], 2 );
              distance = sqrt( distance );
              if ( distance <= radius )
              {
                anyvalid = 1;
                flagbuf[i * numitemsino + j] = 1;
                distancebuf[i * numitemsino + j] = distance;
              }
            }
          }
          else
          {
            anyvalid = 1;
            flagbuf[i * numitemsino + j] = 1;
            distancebuf[i * numitemsino + j] = distance;
          }
        }
      }
      if ( !anyvalid )
      {
        if ( out_invalidhandle )
          DXSetElementInvalid( out_invalidhandle, i );
        *anyinvalid = 1;
      }
    }
  }

  else if ( shape_base == 3 )
  {
    for ( i = 0; i < numitemsbase; i++ )
    {
      anyvalid = 0;
      for ( j = 0; j < numitemsino; j++ )
      {
        if ( DXIsElementValid( in_invalidhandle, j ) )
        {
          if ( radius != -1 )
          {
            if ( ( ( pos_ino_ptr[3 * j] - pos_base_ptr[3 * i] ) < radius ) &&
                 ( ( pos_ino_ptr[3 * j + 1] - pos_base_ptr[3 * i + 1] ) <
                   radius ) &&
                 ( ( pos_ino_ptr[3 * j + 2] - pos_base_ptr[3 * i + 2] ) <
                   radius ) )
            {
              distance =
                  POW( pos_ino_ptr[3 * j] - pos_base_ptr[3 * i], 2 ) +
                  POW( pos_ino_ptr[3 * j + 1] - pos_base_ptr[3 * i + 1], 2 ) +
                  POW( pos_ino_ptr[3 * j + 2] - pos_base_ptr[3 * i + 2], 2 );
              distance = sqrt( distance );
              if ( distance <= radius )
              {
                anyvalid = 1;
                flagbuf[i * numitemsino + j] = 1;
                distancebuf[i * numitemsino + j] = distance;
              }
            }
          }
          else
          {
            anyvalid = 1;
            flagbuf[i * numitemsino + j] = 1;
            distancebuf[i * numitemsino + j] = distance;
          }
        }
      }
      if ( !anyvalid )
      {
        if ( out_invalidhandle )
          DXSetElementInvalid( out_invalidhandle, 1 );
        *anyinvalid = 1;
      }
    }
  }
  else
  {
    DXSetError( ERROR_BAD_PARAMETER, "only 1D 2D 3D positions accepted" );
    goto error;
  }
  return OK;
error:
  return ERROR;
}

/*Added by Jeff Braun to put scattered data on closest grid point*/

Error _dxfConnectScatterObject( Object ino, Object base, Array missing )
{
  struct arg_scatter arg;
  Object subbase;
  int i;

  switch ( DXGetObjectClass( base ) )
  {
    case CLASS_FIELD:
      /* create the task group */
      arg.ino = ino;
      arg.base = (Field)base;
      arg.missing = missing;
      if ( !DXAddTask( ConnectScatterField, ( Pointer ) & arg, sizeof( arg ),
                       0.0 ) )
        goto error;
      break;
    case CLASS_GROUP:
      for ( i = 0; ( subbase = DXGetEnumeratedMember( (Group)base, i, NULL ) );
            i++ )
      {
        if ( !( _dxfConnectScatterObject( (Object)ino, (Object)subbase,
                                          missing ) ) )
          goto error;
      }
      break;
    default:
      DXSetError( ERROR_DATA_INVALID, "base must be a field or a group" );
      goto error;
  }
  return OK;

error:
  return ERROR;
}

static Error ConnectScatterField( Pointer p )
{
  struct arg_scatter *arg = (struct arg_scatter *)p;
  Object ino, newino = NULL;
  Field base;
  Array missing;

  ino = arg->ino;
  newino = DXApplyTransform( ino, NULL );

  switch ( DXGetObjectClass( newino ) )
  {

    case CLASS_FIELD:
      base = arg->base;
      missing = arg->missing;
      if ( !ConnectScatter( (Field)newino, base, missing ) )
        goto error;
      break;

    case CLASS_GROUP:
      DXSetError( ERROR_NOT_IMPLEMENTED, "input cannot be a group" );
      goto error;
      break;

    default:
      DXSetError( ERROR_DATA_INVALID, "input must be a field" );
      goto error;
  }
  DXDelete( (Object)newino );
  return OK;

error:
  DXDelete( (Object)newino );
  return ERROR;
}

static Error ConnectScatter( Field ino, Field base, Array missing )
{
  Matrix matrix, inverse;
  Matrix2D matrix2D, inverse2D;
  Type type;
  Category category;
  Array positions_base, positions_ino, oldcomponent = NULL, newcomponent;
  InvalidComponentHandle out_invalidhandle = NULL, in_invalidhandle = NULL;

  int rank, shape[8];
  int counts[8], rank_base, shape_base[8], shape_ino[8], rank_ino;
  int nitemsbase, nitemsino, nitems, compcount, componentsdone;
  float *ino_positions;
  char *name;

  int i, j, k, ijk[3], index, *count = NULL;
  float xyz[3], sum, eps = .000001;

  float *old_data_f = NULL, *new_data_ptr_f = NULL, *missingval_f = NULL;
  int *old_data_i = NULL, *new_data_ptr_i = NULL, *missingval_i = NULL;
  short *old_data_s = NULL, *new_data_ptr_s = NULL, *missingval_s = NULL;
  double *old_data_d = NULL, *new_data_ptr_d = NULL, *missingval_d = NULL;
  byte *old_data_b = NULL, *new_data_ptr_b = NULL, *missingval_b = NULL;
  ubyte *old_data_ub = NULL, *new_data_ptr_ub = NULL, *missingval_ub = NULL;
  ushort *old_data_us = NULL, *new_data_ptr_us = NULL, *missingval_us = NULL;
  uint *old_data_ui = NULL, *new_data_ptr_ui = NULL, *missingval_ui = NULL;

  if ( DXEmptyField( (Field)ino ) )
    return OK;

  positions_base = (Array)DXGetComponentValue( (Field)base, "positions" );
  if ( !positions_base )
  {
    DXSetError( ERROR_MISSING_DATA, "base field has no positions" );
    goto error;
  }
  positions_ino = (Array)DXGetComponentValue( (Field)ino, "positions" );
  if ( !positions_ino )
  {
    DXSetError( ERROR_MISSING_DATA, "input field has no positions" );
    goto error;
  }

  /*Grid Position Data*/
  DXGetArrayInfo( positions_base, &nitemsbase, &type, &category, &rank_base,
                  shape_base );
  if ( ( DXQueryGridPositions( (Array)positions_base, NULL, counts,
                               (float *)&matrix.b, (float *)&matrix.A ) ) )
  {
    if ( shape_base[0] == 2 )
    {
      matrix2D.A[0][0] = matrix.A[0][0];
      matrix2D.A[1][0] = matrix.A[0][1];
      matrix2D.A[0][1] = matrix.A[0][2];
      matrix2D.A[1][1] = matrix.A[1][0];
      matrix2D.b[0] = matrix.b[0];
      matrix2D.b[1] = matrix.b[1];
      inverse2D = Invert2D( matrix2D );
      inverse.A[0][0] = inverse2D.A[0][0];
      inverse.A[0][1] = inverse2D.A[0][1];
      inverse.A[1][0] = inverse2D.A[1][0];
      inverse.A[1][1] = inverse2D.A[1][1];
    }
    else if ( shape_base[0] == 1 )
    {
      /*origin = matrix.b[0];*/
      /*delta = matrix.A[0][0];*/
      DXSetError( ERROR_DATA_INVALID, "Only 2D and 3D grids supported" );
      goto error;
    }
  }
  else
  {
    /* irregular positions */
    DXSetError( ERROR_DATA_INVALID,
                "irregular grids not currently supported for radius=0" );
    goto error;
  }

  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "positions component must be type float, category real" );
    goto error;
  }

  /*Scatter Position Data*/
  DXGetArrayInfo( positions_ino, &nitemsino, &type, &category, &rank_ino,
                  shape_ino );

  if ( ( type != TYPE_FLOAT ) || ( category != CATEGORY_REAL ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "positions component must be type float, category real" );
    goto error;
  }

  if ( rank_base == 0 )
    shape_base[0] = 1;
  if ( rank_base > 1 )
  {
    DXSetError( ERROR_DATA_INVALID,
                "rank of base positions cannot be greater than 1" );
    goto error;
  }
  if ( rank_ino == 0 )
    shape_ino[0] = 1;
  if ( rank_ino > 1 )
  {
    DXSetError( ERROR_DATA_INVALID,
                "rank of input positions cannot be greater than 1" );
    goto error;
  }

  if ( ( shape_base[0] < 0 ) || ( shape_base[0] > 3 ) )
  {
    DXSetError( ERROR_DATA_INVALID,
                "only 1D, 2D, and 3D positions for base supported" );
    goto error;
  }

  if ( shape_base[0] != shape_ino[0] )
  {
    DXSetError( ERROR_DATA_INVALID, "dimensionality of base (%dD) does not "
                                    "match dimensionality of input (%dD)",
                shape_base[0], shape_ino[0] );
    goto error;
  }

  /* first set the extra counts to one to simplify the code for handling
     1, 2, and 3 dimensions */

  if ( shape_base[0] == 1 )
  {
    counts[1] = 1;
    counts[2] = 1;
  }
  else if ( shape_base[0] == 2 )
  {
    counts[2] = 1;
  }

  /*Will mark as invalid if 'missing' not specified*/

  if ( !missing )
  {
    out_invalidhandle =
        DXCreateInvalidComponentHandle( (Object)base, NULL, "positions" );
    if ( !out_invalidhandle )
      goto error;
    DXSetAllInvalid( out_invalidhandle );
  }
  in_invalidhandle =
      DXCreateInvalidComponentHandle( (Object)ino, NULL, "positions" );
  if ( !in_invalidhandle )
    return ERROR;

  /* need to get the components we are going to copy into the new field */
  /* copied portions from ConnectNearest -JAB */
  compcount = 0;
  componentsdone = 0;
  while ( NULL != ( oldcomponent = (Array)DXGetEnumeratedComponentValue(
                        ino, compcount, &name ) ) )
  {
    Object attr;

    if ( DXGetObjectClass( (Object)oldcomponent ) != CLASS_ARRAY )
      goto component_done;

    if ( !strcmp( name, "positions" ) || !strcmp( name, "connections" ) ||
         !strcmp( name, "invalid positions" ) )
      goto component_done;

    /* if the component refs anything, leave it */
    if ( DXGetComponentAttribute( (Field)ino, name, "ref" ) )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "der" );
    if ( attr )
      goto component_done;

    attr = DXGetComponentAttribute( (Field)ino, name, "dep" );
    if ( !attr )
    {
      /* does it der? if so, throw it away */
      attr = DXGetComponentAttribute( (Field)ino, name, "der" );
      if ( attr )
        goto component_done;
      /* it doesn't ref, doesn't dep and doesn't der.
         Copy it to the output and quit */
      if ( !DXSetComponentValue( (Field)base, name, (Object)oldcomponent ) )
        goto error;
      goto component_done;
    }

    if ( DXGetObjectClass( attr ) != CLASS_STRING )
    {
      DXSetError( ERROR_DATA_INVALID, "dependency attribute not type string" );
      goto error;
    }

    if ( strcmp( DXGetString( (String)attr ), "positions" ) )
    {
      DXWarning( "Component '%s' is not dependent on positions!", name );
      DXWarning( "Regrid will remove '%s' and output the base grid 'data' "
                 "component if it exists",
                 name );
      goto component_done;
    }

    /* if the component happens to be "invalid positions" ignore it */
    if ( !strcmp( name, "invalid positions" ) )
      goto component_done;

    componentsdone++;

    DXGetArrayInfo( (Array)oldcomponent, &nitems, &type, &category, &rank,
                    shape );
    if ( rank == 0 )
      shape[0] = 1;

    /* check that the missing component, if present, matches the component
       we're working on. If missing not present, 0 is used as default */
    switch ( type )
    {
      case( TYPE_FLOAT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    float, type, shape[0], missingval_f );
        break;
      case( TYPE_INT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ), int,
                                    type, shape[0], missingval_i );
        break;
      case( TYPE_DOUBLE ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    double, type, shape[0], missingval_d );
        break;
      case( TYPE_SHORT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    short, type, shape[0], missingval_s );
        break;
      case( TYPE_BYTE ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    byte, type, shape[0], missingval_b );
        break;
      case( TYPE_UINT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    uint, type, shape[0], missingval_ui );
        break;
      case( TYPE_USHORT ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    ushort, type, shape[0], missingval_us );
        break;
      case( TYPE_UBYTE ) :
        CHECK_AND_ALLOCATE_MISSING( missing, DXGetItemSize( oldcomponent ),
                                    ubyte, type, shape[0], missingval_ub );
        break;
      default:
        DXSetError( ERROR_DATA_INVALID, "unsupported data type" );
        goto error;
    }

    if ( rank > 1 )
    {
      DXSetError( ERROR_NOT_IMPLEMENTED, "rank larger than 1 not implemented" );
      goto error;
    }
    if ( nitems != nitemsino )
    {
      DXSetError(
          ERROR_BAD_PARAMETER,
          "number of items in %s component not equal to number of positions",
          name );
      goto error;
    }

    /*Intialize newcomponent with missing value*/
    newcomponent = DXNewArrayV( type, category, rank, shape );
    if ( !DXAddArrayData( newcomponent, 0, nitemsbase, NULL ) )
      goto error;

    switch ( type )
    {
      case( TYPE_FLOAT ) :
        old_data_f = (float *)DXGetArrayData( oldcomponent );
        new_data_ptr_f = (float *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_f[i * shape[0] + j] = missingval_f[j];
          }
        }
        break;
      case( TYPE_INT ) :
        old_data_i = (int *)DXGetArrayData( oldcomponent );
        new_data_ptr_i = (int *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_i[i * shape[0] + j] = missingval_i[j];
          }
        }
        break;
      case( TYPE_SHORT ) :
        old_data_s = (short *)DXGetArrayData( oldcomponent );
        new_data_ptr_s = (short *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_s[i * shape[0] + j] = missingval_s[j];
          }
        }
        break;
      case( TYPE_BYTE ) :
        old_data_b = (byte *)DXGetArrayData( oldcomponent );
        new_data_ptr_b = (byte *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_b[i * shape[0] + j] = missingval_b[j];
          }
        }
        break;
      case( TYPE_DOUBLE ) :
        old_data_d = (double *)DXGetArrayData( oldcomponent );
        new_data_ptr_d = (double *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_d[i * shape[0] + j] = missingval_d[j];
          }
        }
        break;
      case( TYPE_UINT ) :
        old_data_ui = (uint *)DXGetArrayData( oldcomponent );
        new_data_ptr_ui = (uint *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_ui[i * shape[0] + j] = missingval_ui[j];
          }
        }
        break;
      case( TYPE_UBYTE ) :
        old_data_ub = (ubyte *)DXGetArrayData( oldcomponent );
        new_data_ptr_ub = (ubyte *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_ub[i * shape[0] + j] = missingval_ub[j];
          }
        }
        break;
      case( TYPE_USHORT ) :
        old_data_us = (ushort *)DXGetArrayData( oldcomponent );
        new_data_ptr_us = (ushort *)DXGetArrayData( newcomponent );
        for ( i = 0; i < nitemsbase; i++ )
        {
          for ( j = 0; j < shape[0]; j++ )
          {
            new_data_ptr_us[i * shape[0] + j] = missingval_us[j];
          }
        }
        break;
      default:
        DXSetError( ERROR_DATA_INVALID, "unrecognized data type" );
        goto error;
    }

    /*Need to allocate count and compute inverse of grid matrix*/
    count = (int *)DXAllocateLocalZero( sizeof( int ) * nitemsbase );
    if ( !count )
      goto error;
    if ( shape_base[0] == 3 )
    {
      inverse = DXInvert( matrix );
    }

    ino_positions = (float *)DXGetArrayData( positions_ino );

    /*Put Scattered data into grid data array*/
    for ( i = 0; i < nitemsino; i++ )
    {/*loops through scatter data */
      if ( DXIsElementValid( in_invalidhandle, i ) )
      {/*seems to work fine without this*/
        index = 0;
        for ( j = 0; j < shape_base[0]; j++ )
        {
          xyz[j] = ino_positions[shape_base[0] * i + j] - matrix.b[j];
        }
        for ( j = 0; j < shape_base[0]; j++ )
        {
          sum = 0;
          for ( k = 0; k < shape_base[0]; k++ )
          {
            sum = sum + inverse.A[k][j] * xyz[k];
          }
          ijk[j] = rint( sum + eps );
          if ( ijk[j] < 0 || ijk[j] >= counts[j] )
          {
            index = -1;
          }
        }

        /* Finds data array index number for grid position*/
        if ( index != -1 )
        {
          for ( j = 0; j < shape_base[0]; j++ )
          {
            for ( k = shape_base[0] - 1; k > j; k-- )
            {
              ijk[j] = ijk[j] * counts[k];
            }
            index = index + ijk[j];
          }
        }

        if ( index >= 0 && index < nitemsbase )
        {
          if ( !missing )
          {
            if ( !DXSetElementValid( out_invalidhandle, index ) )
              goto error;
          }
          switch ( type )
          {
            case( TYPE_FLOAT ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_f[shape[0] * index + j] =
                      old_data_f[shape[0] * i + j];
                else
                  new_data_ptr_f[shape[0] * index + j] =
                      ( new_data_ptr_f[shape[0] * index + j] * count[index] +
                        old_data_f[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_INT ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_i[shape[0] * index + j] =
                      old_data_i[shape[0] * i + j];
                else
                  new_data_ptr_i[shape[0] * index + j] =
                      ( new_data_ptr_i[shape[0] * index + j] * count[index] +
                        old_data_i[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_SHORT ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_s[shape[0] * index + j] =
                      old_data_s[shape[0] * i + j];
                else
                  new_data_ptr_s[shape[0] * index + j] =
                      ( new_data_ptr_s[shape[0] * index + j] * count[index] +
                        old_data_s[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_DOUBLE ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_d[shape[0] * index + j] =
                      old_data_d[shape[0] * i + j];
                else
                  new_data_ptr_d[shape[0] * index + j] =
                      ( new_data_ptr_d[shape[0] * index + j] * count[index] +
                        old_data_d[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_BYTE ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_b[shape[0] * index + j] =
                      old_data_b[shape[0] * i + j];
                else
                  new_data_ptr_b[shape[0] * index + j] =
                      ( new_data_ptr_b[shape[0] * index + j] * count[index] +
                        old_data_b[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_UINT ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_ui[shape[0] * index + j] =
                      old_data_ui[shape[0] * i + j];
                else
                  new_data_ptr_ui[shape[0] * index + j] =
                      ( new_data_ptr_ui[shape[0] * index + j] * count[index] +
                        old_data_ui[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_USHORT ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_us[shape[0] * index + j] =
                      old_data_us[shape[0] * i + j];
                else
                  new_data_ptr_us[shape[0] * index + j] =
                      ( new_data_ptr_us[shape[0] * index + j] * count[index] +
                        old_data_us[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            case( TYPE_UBYTE ) :
              for ( j = 0; j < shape[0]; j++ )
              {
                if ( count[index] == 0 )
                  new_data_ptr_ub[shape[0] * index + j] =
                      old_data_ub[shape[0] * i + j];
                else
                  new_data_ptr_ub[shape[0] * index + j] =
                      ( new_data_ptr_ub[shape[0] * index + j] * count[index] +
                        old_data_ub[shape[0] * i + j] ) /
                      ( count[index] + 1 );
              }
              break;
            default:
              DXSetError( ERROR_DATA_INVALID, "unrecognized data type" );
              goto error;
          }
          count[index]++;
        }
      }
    }
    if ( !DXChangedComponentValues( base, name ) )
      goto error;
    if ( !DXSetComponentValue( base, name, (Object)newcomponent ) )
      goto error;
    newcomponent = NULL;
    if ( !DXSetComponentAttribute( (Field)base, name, "dep",
                                   (Object)DXNewString( "positions" ) ) )
      goto error;

  component_done:
    compcount++;
    DXFree( (Pointer)missingval_f );
    missingval_f = NULL;
    DXFree( (Pointer)missingval_s );
    missingval_s = NULL;
    DXFree( (Pointer)missingval_i );
    missingval_i = NULL;
    DXFree( (Pointer)missingval_d );
    missingval_d = NULL;
    DXFree( (Pointer)missingval_b );
    missingval_b = NULL;
    DXFree( (Pointer)missingval_ui );
    missingval_ui = NULL;
    DXFree( (Pointer)missingval_ub );
    missingval_ub = NULL;
    DXFree( (Pointer)missingval_us );
    missingval_us = NULL;
    DXFree( (Pointer)count );
    count = NULL;
  }

  if ( !missing )
  {
    for ( i = 0; i < nitemsbase; i++ )
    {
      if ( DXIsElementInvalid( out_invalidhandle, i ) )
      {
        i = nitemsbase;
        if ( !DXSaveInvalidComponent( base, out_invalidhandle ) )
          goto error;
      }
    }
  }
  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  if ( !DXEndField( base ) )
    goto error;
  return OK;

error:
  DXFreeInvalidComponentHandle( out_invalidhandle );
  DXFreeInvalidComponentHandle( in_invalidhandle );
  DXFree( (Pointer)missingval_f );
  DXFree( (Pointer)missingval_s );
  DXFree( (Pointer)missingval_i );
  DXFree( (Pointer)missingval_d );
  DXFree( (Pointer)missingval_b );
  DXFree( (Pointer)missingval_ui );
  DXFree( (Pointer)missingval_ub );
  DXFree( (Pointer)missingval_us );
  DXFree( (Pointer)count );
  return ERROR;
}
